#include "pvmodel_GPU.cuh"

#include <cuda_runtime.h>
#include <cublas_v2.h>
#include "device_launch_parameters.h"
#include "cusolverDn.h"

#include <ctime>

__global__ void Control_Init(double* time, int* Dblk, int* SW_Char)
{
    int idx = threadIdx.x;

    if(*time > 0.08)
    {
        Dblk[idx] = 1;
    }

    if(*time > 0.05)
    {  
        SW_Char[idx] = 1;
    }
}

__device__ void LPF_Execute(double Kf, double G, double input, double* output_his, double* input_his, double* output)
{
    *output = 1.0 / (1 + Kf) * G * (*input_his + input) + (Kf - 1) / (Kf + 1) * (*output_his);
    *output_his = *output;
    *input_his = input;
}

__device__ void signal_switch_Execute(double input1, double input2, int ctrl, int signal, double* output)
{
    if (signal == ctrl)
        *output = input1;  
    else
        *output = input2;
}

__device__ void rate_limiter_Execute(double input, double dt, double ULIM, double LLIM, double* output_his, double* output)
{   
    double rate = (input - *output_his) / dt;
    if (rate > 0)
    {
        if (rate > ULIM)
        {
            *output = *output_his + ULIM * dt;
        }
        else
        {
            *output = input;
        }
    }
    else
    {
        if (fabs(rate) >= LLIM)
        {
            *output = *output_his - LLIM * dt;
        }
        if (fabs(rate) < LLIM)
        {
            *output = input;
        }
    }
    
    *output_his = *output;
}

__device__ void PI_model_Execute(PI_model* PImodel, double PI_IN, int RST, int FRZ, double MAXOUT, double MINOUT, double T, double Kp, double DELT)
{
    PImodel->T1 = DELT / (2 * T);
    
    //Not Gate
    if (RST == 1)
        PImodel->k1 = 0;
    else
        PImodel->k1 = 1;

    //And Gate
    if (FRZ == 1 || PImodel->Y5_his == 1 || PImodel->Y4_his == 1)
        PImodel->k4 = 1;
    else
        PImodel->k4 = 0;

    //Comparator
    if (PImodel->Y9_his > MAXOUT)
        PImodel->k3 = 1;
    else
        PImodel->k3 = 0;

    //Comparator
    if (MINOUT > PImodel->Y9_his)
        PImodel->k2 = 1;
    else
        PImodel->k2 = 0;

    //signal switch
    if (PImodel->Y6_his == 1)
        PImodel->k5 = 0;
    else
        PImodel->k5 = 1;

    //Limiter
    if (PImodel->Y9_his > MAXOUT)
    {
        PImodel->k6 = 1;
        PImodel->k7 = PImodel->k8 = 0;
    }
    else if ((PImodel->Y9_his < MINOUT))
    {
        PImodel->k6 = PImodel->k8 = 0;
        PImodel->k7 = 1;
    }
    else
    {
        PImodel->k8 = 1;
        PImodel->k6 = PImodel->k7 = 0;
    }

    //Integrator Limiter
    if (PImodel->Y8_his > MAXOUT)
    {
        PImodel->k9 = 1;
        PImodel->k10 = PImodel->k11 = 0;
    }
    else if ((PImodel->Y8_his < MINOUT))
    {
        PImodel->k10 = 1;
        PImodel->k9 = PImodel->k11 = 0;
    }
    else
    {
        PImodel->k11 = 1;
        PImodel->k9 = PImodel->k10 = 0;
    }

    PImodel->Y1 = PImodel->k1;
    PImodel->Y2 = PI_IN * PImodel->Y1;
    PImodel->Y3 = Kp * PImodel->Y2;
    PImodel->Y4 = PImodel->k2;
    PImodel->Y5 = PImodel->k3;
    PImodel->Y6 = PImodel->k4;
    PImodel->Y7 = PImodel->k5 * PI_IN;
    PImodel->Y8 = PImodel->k9 * MAXOUT + PImodel->k10 * MINOUT + PImodel->k11 * (PImodel->Y8_his + PImodel->T1 * (PImodel->Y7 + PImodel->Y7_his));

    PImodel->Y9 = PImodel->Y3 + PImodel->Y8;
    PImodel->Y10 =  PImodel->k8 * Kp *  PImodel->k1 * PI_IN +  PImodel->k8 *  PImodel->Y8 +  PImodel->k6 * MAXOUT +  PImodel->k7 * MINOUT;

    //output
    PImodel->PI_OUT = PImodel->Y10;

    //History Value
    PImodel->Y1_his = PImodel->Y1;
    PImodel->Y4_his = PImodel->Y4;
    PImodel->Y5_his = PImodel->Y5;
    PImodel->Y6_his = PImodel->Y6;

    PImodel->Y7_his = PImodel->Y7;
    PImodel->Y8_his = PImodel->Y8;
    PImodel->Y9_his = PImodel->Y9;

}

__device__ void rec_polar_Execute(double input_X, double input_Y, double* output_M, double* output_P)
{
    *output_M = sqrt(pow(input_X, 2) + pow(input_Y, 2));
    *output_P = atan2(input_Y, input_X);
}

__device__ void limiter_Execute(double input, double MAX, double MIN, double * output)
{
    if (input > MAX)
    {
        *output = MAX;
    }
    else if (input < MIN)
    {
        *output = MIN;
    }
    else
    {
        *output = input;
    }
}

__device__ void polar_rec_Execute(double input_M, double input_P, double* output_X, double* output_Y)
{
    *output_X = input_M * cos(input_P);
    *output_Y = input_M * sin(input_P);
}

__device__ void Two_PI_Limit_Execute(double W, double* PHI1)
{
    double W1 = W;
    if (W1 < 2.0 * PI)
    {
        if (W1 > 0)
        {
            *PHI1 = W1;
        }
        else
        {
            while (W1 < 0)
            {
                W1 = W1 + 2.0 * PI;
            }
            *PHI1 = W1;
        }
    }
    else
    {
        while (W1 > 2.0 * PI)
        {
            W1 = W1 - 2.0 * PI;
        }
        *PHI1 = W1;
    }

}


__device__ void PLL3Phase_Execute(PLL3Phase* PLL, double KP, double KI, double BASE, double BASEFREQ, double fmax, double fmin, double deltaT, double inputA, double inputB, double inputC)
{
    //Init
    PLL->PIPG = KP;
    PLL->PIIT = 1.0 / KI;
    PLL->Base = BASE;
    PLL->Frat = BASEFREQ;
    PLL->dt = deltaT;
    PLL->fc = 70.5;
    PLL->Fmax = fmax * PLL->Frat;
    PLL->Fmin = fmin * PLL->Frat;
    
    PLL->Fmax1 = PLL->Fmax * 2.0 * PI;
    PLL->Fmin1 = PLL->Fmin * 2.0 * PI;
    
    PLL->Tf = 1.0 / (2.0 * PI * PLL->fc);
    PLL->Kf = 2 * PLL->Tf / PLL->dt;

    // Clark_transformation
    double a, b;
    a = 2.0 / 3 * (inputA - inputB / 2 - inputC / 2) / PLL->Base;
    b = sqrt(3.0) / 3 * (inputB - inputC) / PLL->Base;

    //Park_transformation
    double vd1, vq1;
    vd1 = a * cos(PLL->PHI_his) + b * sin(PLL->PHI_his);
    vq1 = -a * sin(PLL->PHI_his) + b * cos(PLL->PHI_his);

    //iPark_transformation
    double vd2, vq2;
    vd2 = a * cos(PLL->PHI_his) - b * sin(PLL->PHI_his);
    vq2 = a * sin(PLL->PHI_his) + b * cos(PLL->PHI_his);

    double sumV = inputA + inputB + inputC;

    double Err_integral;
    double PI_out;

    if (sumV <= 1e-8 && sumV >= -1e-8) // three-phase equilibrium//
    {
        // PI_Backward_Euler//
        Err_integral = (PLL->dt / PLL->PIIT * 0.5) * (vq1 + PLL->Err_input_old) + PLL->Err_integral_old;
        PLL->Err_input_old = vq1;
        PLL->Err_integral_old = Err_integral;

        PI_out = PLL->PIPG * vq1 + Err_integral + 2.0 * PI * PLL->Frat;
    }
    else // three-phase imbalance//
    {

        double output_LPF1, output_LPF2, output_LPF3, output_LPF4;
        LPF_Execute(PLL->Kf, 1.0, PLL->vd1_1_his, &(PLL->LPF1_output_his), &(PLL->LPF1_input_his), &output_LPF1);
        LPF_Execute(PLL->Kf, 1.0, PLL->vq1_1_his, &(PLL->LPF2_output_his), &(PLL->LPF2_input_his), &output_LPF2);
        LPF_Execute(PLL->Kf, 1.0, PLL->vd2_1_his, &(PLL->LPF3_output_his), &(PLL->LPF3_input_his), &output_LPF3);
        LPF_Execute(PLL->Kf, 1.0, PLL->vq2_1_his, &(PLL->LPF4_output_his), &(PLL->LPF4_input_his), &output_LPF4);

        double vd1_L = output_LPF1;
        double vq1_L = output_LPF2;
        double vd2_L = output_LPF3;
        double vq2_L = output_LPF4;

        double vd1_1 = vd1 - vd2_L * cos(2 * PLL->PHI_his) - vq2_L * sin(2 * PLL->PHI_his);
        double vq1_1 = vq1 + vd2_L * sin(2 * PLL->PHI_his) - vq2_L * cos(2 * PLL->PHI_his);
        double vd2_1 = vd2 - vd1_L * cos(2 * PLL->PHI_his) + vq1_L * sin(2 * PLL->PHI_his);
        double vq2_1 = vq2 - vd1_L * sin(2 * PLL->PHI_his) - vq1_L * cos(2 * PLL->PHI_his);

        PLL->vd1_1_his = vd1_1;
        PLL->vq1_1_his = vq1_1;
        PLL->vd2_1_his = vd2_1;
        PLL->vq2_1_his = vq2_1;

        // PI_Backward_Euler
        Err_integral = (PLL->dt / PLL->PIIT * 0.5) * (vq1_1 + PLL->Err_input_old) + PLL->Err_integral_old;
        PLL->Err_input_old = vq1_1;
        PLL->Err_integral_old = Err_integral;

        PI_out = PLL->PIPG * vq1_1 + Err_integral + 2.0 * PI * PLL->Frat;
    }

    double W1;
    limiter_Execute(PI_out, PLL->Fmax1, PLL->Fmin1, &W1);

    //integral
    double output_integral = (PLL->input_integral_old + W1) * (PLL->dt / 2.0) + PLL->output_integral_old;
    PLL->input_integral_old = W1;
    PLL->output_integral_old = output_integral;

    output_integral = output_integral - 2.0 * PI;

    // 0_2PI_Limiter
    double output_y;
    Two_PI_Limit_Execute(output_integral, &output_y);

    PLL->output = output_y;
    PLL->PHI_his = PLL->output;

}

__device__ void angle_resolver_Execute(double input, double* output)
{
    double n;
        
    if (input > 2.0 * PI)
    {
        n = int(input / (2.0 * PI));
        *output = input - 2.0 * PI * n;
    }
    else if (input < 0)
    {
        n = int(fabs(input) / (2.0 * PI));
        *output = input + (n + 1) * 2.0 * PI;
    }
    else
    {
        *output = input;
    }
}

__device__ void ABC_DQO_Execute(double abc0, double abc1, double abc2, double theta, double* output0, double* output1)
{
    double dq00 = (2.0 / 3.0) * (cos(theta) * abc0 + cos(theta - 2.0 * acos(-1.0) / 3.0) * abc1 + cos(theta + 2.0 * acos(-1.0) / 3.0) * abc2);
    double dq01 = (2.0 / 3.0) * (sin(theta) * abc0 + sin(theta - 2.0 * acos(-1.0) / 3.0) * abc1 + sin(theta + 2.0 * acos(-1.0) / 3.0) * abc2);
    *output0 = dq00;
    *output1 = dq01;
}

__device__ void DQ0_ABC_Execute(double dq0, double dq1, double dq2, double theta, double* output0, double* output1, double* output2)
{
    double abc0 = cos(theta) * dq1 - sin(theta) * dq2 + dq0;
    double abc1 = cos(theta - 2.0 * acos(-1.0) / 3.0) * dq1 - sin(theta - 2.0 * acos(-1.0) / 3) * dq2 + dq0;
    double abc2 = cos(theta + 2.0 * acos(-1.0) / 3.0) * dq1 - sin(theta + 2.0 * acos(-1.0) / 3.0) * dq2 + dq0;
    *output0 = abc0;
    *output1 = abc1;
    *output2 = abc2;
}

__device__ void comparator_Execute(double input_A, double input_B, double comparator_H, double comparator_L, double* output)
{
    if (input_A > input_B)
    {
        *output = comparator_H;
    }
    else
    {
        *output = comparator_L;
    }
}

__device__ void SR_Flip_Flop_Execute(int input_s, int input_r, int* input_s_his, int* input_r_his, int* output_Q_his, int* output_Q_not_his, int* output_Q, int* output_Q_not)
{
    if ((input_s > *input_s_his) || (input_r > *input_r_his))
    {
        if (input_s == 0)
        {
            *output_Q_not = 0;
            *output_Q = 1;
        }
        else
        {
            *output_Q_not = 1;
            *output_Q = 0;
        }
           
    }
    else
    {
        *output_Q = *output_Q_his;
        *output_Q_not = *output_Q_not_his;
    }

    if ((input_s < *input_s_his) || (input_r < *input_r_his))
    {
            
        if (input_s == 0)
        {
            *output_Q_not = 0;
            *output_Q = 1;
        }
        else
        {
            *output_Q_not = 1;
            *output_Q = 0;
        }
           
    }
    else
    {
        *output_Q = *output_Q_his;
        *output_Q_not = *output_Q_not_his;
    }

    *input_r_his = input_r;
    *input_s_his = input_s;
    *output_Q_his = *output_Q;
    *output_Q_not_his = *output_Q_not;

}

__device__ void Two_input_and_gate_Execute(int input1, int input2, int* output)
{
    int A[2], B[2];
    int i;
    int out;

    A[0] = abs(input1);
    A[1] = abs(input2);

    for (i = 0; i < 2; i++)
    {
        B[i] = A[i] % 2;
        A[i] = B[i];
    }
    out = 1;
    for (i = 0; i < 2; i++)
    {
        out = out & A[i];
    }
    *output = out;
}

__device__ void MPPT_Execute(MPPT* mppt, int Reset, double Vinput, double Iinput, double Sd, double Vstep)
{

    if ((Reset - mppt->Reset_old) == 1)
    {
        mppt->output_State = 0;
    }

    mppt->Reset_old = Reset;

    if (mppt->SampleCount > 0)
    {
        mppt->SampleCount = mppt->SampleCount - 1;
    }

  
    
    if (mppt->SampleCount == 0)
    {
        if (mppt->output_State == 0)
        {
            if (mppt->StartBlk == 0)
            {
                mppt->output_Vopen = Vinput;
                mppt->output_VstartStep = mppt->output_Vopen / mppt->output_StrtupNum;
                mppt->StartBlk = 1;
                mppt->output_Vmppt = 0;
            }

            mppt->output_Vmppt = mppt->output_Vmppt + mppt->output_VstartStep;
            mppt->StrtupCount = mppt->StrtupCount - 1;

            if (mppt->StrtupCount == 0)
            {
                mppt->output_State = 1;
                mppt->StrtupCount = mppt->output_StrtupNum;
                mppt->StartBlk = 0;
            }
        }

        if (mppt->output_State == 1)
        {
            mppt->Vmppt_old = mppt->output_Vmppt;
            double dI = Iinput - mppt->Iold;
            double dV = Vinput - mppt->Vold;
            double dP = Iinput * Vinput - mppt->Iold * mppt->Vold;

            mppt->Iold = Iinput;
            mppt->Vold = Vinput;

            if (dV != 0)
            {
                mppt->Cond = (dI / dV) * Vinput + Iinput;
                if (mppt->Cond > Sd)
                    mppt->output_Vmppt = mppt->output_Vmppt + Vstep * 0.001;

                if (mppt->Cond < (-1.0 * Sd))
                    mppt->output_Vmppt = mppt->output_Vmppt - Vstep * 0.001;
            }
            else
            {
                if (dI > 0.0)
                    mppt->output_Vmppt = mppt->output_Vmppt + Vstep * 0.001;
                
                if (dI < 0.0)
                    mppt->output_Vmppt = mppt->output_Vmppt - Vstep * 0.001;
            }
        }

        if ((mppt->output_Vmppt > mppt->output_Vopen) || (mppt->output_Vmppt < 0.0))
        {
            mppt->output_Vmppt = mppt->Vmppt_old;
        }

        mppt->SampleCount = mppt->output_SampNum;
    
    }
        
        
        //mppt->output = Cond;
        
        mppt->Vmppt_old = mppt->output_Vmppt;
        //Store_Vold = Vold;
        //Store_Iold = Iold;
        //Store_StartBlk = StartBlk;
        //Store_SampleCount = SampleCount;
        //mppt->Reset_old = Reset;
        //Store_StrtupCount = StrtupCount;
}

__device__ void Sawtooth_Generator_Execute(double input_freq, double max_pwm, double min_pwm, double DELT, double* Count_i, double* output)
{
    double NumCount = 1.0 / (input_freq * DELT);
    double slope1 = (max_pwm - (max_pwm + min_pwm) / 2.0) * 4.0;
    double slope2 = -slope1;
    *Count_i = *Count_i + 1.0 / NumCount;
    double slope, intercept;

    while (*Count_i < 0)
    {
        *Count_i = *Count_i + 1;
    }
    
    while (*Count_i > 1)
    {
        *Count_i = *Count_i - 1;
    }

    if (*Count_i >= 0 && *Count_i < 0.25)
    {
        slope = slope1;
        intercept = (max_pwm + min_pwm) / 2.0;
    }
    else if (*Count_i >= 0.25 && *Count_i < 0.75)
    {
        slope = slope2;
        intercept = min_pwm + 3.0 / 2.0 * (max_pwm - min_pwm);
    }
    else
    {
        slope = slope1;
        intercept = min_pwm - (3.0 / 2.0) * (max_pwm - min_pwm);
    }

    *output = slope * (*Count_i) + intercept;
}

__device__ void Interpolated_Firing_Pulses_Execute(int Dblok, double H_OFF, double L_OFF, double H_ON, double L_ON, double delt, double* fx_his, int* output1, double* output2)
{
    if (Dblok == 0)
    {
        *output1 = 0;
        *output2 = 0.0;
    }
    else
    {
        if (H_OFF >= L_OFF)
        {
            *output1 = 0;
        }
        else
        {
            *output1 = 1;
        }
        if (H_ON >= L_ON)
        {
            *output1 = 1;
        }
        else
        {
            *output1 = 0;
        }
        
        double fx = L_OFF - H_OFF;

            //output2 = abs(fx) * delt /(abs(fx) + abs(fx_his));
        if ((*fx_his * fx) <= 0)
        {
            *output2 = fabs(fx) * delt / (fabs(fx) + fabs(*fx_his));
        }
        else
            *output2 = 0.0;


        *fx_his = fx;
    }
    
}

__global__ void Control_PVmodel_kernel(int* Dblk1, double* comparator_1_B, double* comparator_2_A, double* VLLrated, double* fsys, double* Lcon1, double* Vtr2Pri, double* Vtr2Sec, double* PoRated, 
double* Vcon1A, double* Vcon1B, double* Vcon1C, double* Icon1A, double* Icon1B, double* Icon1C, double* VdcBus, double* VdcRated, double* Pcon1filt, double* Qcon1filt, double* VacRMSfilt, double* Ipv, double* VpvBus, 
int* VdcCon, double* VpvRated, double* signal_A, int* FrzPI, double* VdcPURef, double* VacPURef, double* QPURef, int* VacOrQ, double* IdMax, double* IdMin, double* Park0, int* Har3ctrl, double* fsw1, double* fsw2, double* VdcPreCh, 
int* MPPT_ctrl, double* MPPT_ctrl_timetag, int* Ta_1, double* Ta_1_timetag, int* Ta_4, double* Ta_4_timetag, int* Tb_3, double* Tb_3_timetag, int* Tb_6, double* Tb_6_timetag, int* Tc_5, double* Tc_5_timetag, int* Tc_2, double* Tc_2_timetag, 
double* G_VSC, double* T_VSC, double* Strtup, double* Samp, double* Vstep, double* Sd, double* G1_Boost, double* T1_Boost, double* G2_Boost, double* T2_Boost, double* G3_Boost, double* T3_Boost, int* signal, double* Max_Boost, double* Min_Boost, double* Kp_Boost, double* Ti_Boost, double* Kp, double* Kp2, double* ULIM, double* LLIM, double* MAXOUT, double* MINOUT, double* MAXOUT2, double* MINOUT2, double* PI_T1, double* PI_T2, 
double* Mmax_current, double* Mmin_current, double* MaxLimit1_current, double* MinLimit1_current, double* MaxLimit2_current, double* MinLimit2_current, double* Kp1_current, double* Kp2_current, double* Ti1_current, double* Ti2_current, 
double* Kp_PLL, double* Ti_PLL, double* Base_Vpll, double* fbase, double* fmax, double* fmin, double* Kp_realpole, double* Ti_realpole_L, double* Ti_realpole_R, int* Ctrl, double* max1_pwm, double* min1_pwm, double* max2_pwm, double* min2_pwm, 
double* G_start, double* T_start, double* comparator1_H, double* comparator1_L, double* comparator2_H, double* comparator2_L, int* SW_PVunit, int* SW_PV, int* SW_Chopper, double* DELT, double* TIME, double* theta1, double* theta2,
double* state_his, double* output_his_LPF_SELF_start, double* input_his_LPF_SELF_start, int* output_Q_his, int* output_Q_not_his, int* input_s_his, int* input_r_his, double* output_his_LPF_VSC, double* input_his_LPF_VSC,
double* output_his_LPF_PLL, double* input_his_LPF_PLL, PLL3Phase* componentPLL3Phase,
double* output_his_rate_limiter_Voltage, PI_model* componentPI_model_Voltage,
PI_model* componentPI_model_Current,
double* output_his_LPF_boost, double* input_his_LPF_boost, double* output_his_LPF2_boost, double* input_his_LPF2_boost, MPPT* componentMPPT, PI_model* componentPI_model_boost,
double* Count_i_Sawtooth_Generator, double* fx_his_Interpolated_Firing_Pulses)
{
    int idx = threadIdx.x;

    //Start
    int ResetPI2, ResetPI1;
    double StartPI2, StartPI1;
    int Dblk2;

    //CHECK!!!!
    int input_not_gate[3];
    double inputA_comprator[2], inputB_comprator[2];
    double comparator_H[2], comparator_L[2];

    double VdcBusPU;
    VdcBusPU = VdcBus[idx] / *VdcRated;
    
    double output_LPF_SELF_start;
    LPF_Execute(2.0 * (*T_start) / *DELT, *G_start, VdcBusPU, &output_his_LPF_SELF_start[idx], &input_his_LPF_SELF_start[idx], &output_LPF_SELF_start);

    inputA_comprator[0] = inputB_comprator[1] = output_LPF_SELF_start;
    inputA_comprator[1] = *comparator_2_A;
    inputB_comprator[0] = *comparator_1_B;
    comparator_H[0] = *comparator1_H;
    comparator_H[1] = *comparator2_H;
    comparator_L[0] = *comparator1_L;
    comparator_L[1] = *comparator2_H;

    double output_comparator_start[2];
    for(int i = 0; i < 2; i++)
    {
        comparator_Execute(inputA_comprator[i], inputB_comprator[i], comparator_H[i], comparator_L[i], &output_comparator_start[i]);
    }
    
    int output_Q_SR_Flip_Flop, output_Q_not_SR_Flip_Flop;
    SR_Flip_Flop_Execute(output_comparator_start[0], output_comparator_start[1], &input_r_his[idx], &input_s_his[idx], &output_Q_his[idx], &output_Q_not_his[idx], &output_Q_SR_Flip_Flop, &output_Q_not_SR_Flip_Flop);

    Two_input_and_gate_Execute(Dblk1[idx], state_his[idx], &Dblk2);

    input_not_gate[0] = Dblk1[idx];
    input_not_gate[1] = Dblk2;
    input_not_gate[2] = output_Q_SR_Flip_Flop;

    SW_PVunit[idx] = ResetPI1 = 1 - abs(input_not_gate[0]) % 2;
    SW_PV[idx] = ResetPI2 = 1 - abs(input_not_gate[1]) % 2;
    //CHECK!!!!
    SW_Chopper[idx] = 1 - abs(input_not_gate[2]) % 2;

    StartPI1 = double(Dblk1[idx]);
    StartPI2 = double(Dblk2);
    
    //Bases
    double Zcon1, Ntr2;
    double Vtr2PriPh, Vtr2SecPh, VLLcon1, VLLcon1PK, Vphcon1PK;
    double Icon1basePK, wLeqPU;

    Vtr2PriPh = 0.577350 * (*Vtr2Pri);
    Vtr2SecPh = 1.0 * (*Vtr2Sec);
    Ntr2 = Vtr2PriPh / (Vtr2SecPh);

    VLLcon1 = *VLLrated / Ntr2;
    Zcon1 = pow(VLLcon1, 2.0) / (*PoRated);
    VLLcon1PK = 1.41421356 * (VLLcon1);
    Vphcon1PK = 0.57735 * VLLcon1PK;
    VdcPreCh[idx] = 1.1 * VLLcon1PK;
    Icon1basePK = Vphcon1PK / Zcon1;

    wLeqPU = *Lcon1 * (*fsys) * 6.2831853 / Zcon1;

    //VSC
    double divider_input1[10], divider_input2[10];
    double Vcon1ABpu, Vcon1BCpu, Vcon1CApu;
    double Icon1Apu, Icon1Bpu, Icon1Cpu;
    double Pcon1PUfilt, Qcon1PUfilt, VacPUfilt;
    //double VdcBusPU;
    double Kf;

    divider_input1[0] = Vcon1A[idx] - Vcon1B[idx];
    divider_input1[1] = Vcon1B[idx] - Vcon1C[idx];
    divider_input1[2] = Vcon1C[idx] - Vcon1A[idx];
    divider_input2[0] = divider_input2[1] = divider_input2[2] = VLLcon1PK;
    divider_input1[3] = Icon1A[idx];
    divider_input1[4] = Icon1B[idx];
    divider_input1[5] = Icon1C[idx];
    divider_input2[3] = divider_input2[4] = divider_input2[5] = Icon1basePK;
    divider_input1[6] = VdcBus[idx];
    divider_input2[6] = *VdcRated;
    divider_input1[7] = Pcon1filt[idx];
    divider_input2[7] = *PoRated;
    divider_input1[8] = Qcon1filt[idx];
    divider_input2[8] = *PoRated;

    Vcon1ABpu = divider_input1[0] / divider_input2[0];
    Vcon1BCpu = divider_input1[1] / divider_input2[1];
    Vcon1CApu = divider_input1[2] / divider_input2[2];
    Icon1Apu = divider_input1[3] / divider_input2[3];
    Icon1Bpu = divider_input1[4] / divider_input2[4];
    Icon1Cpu = divider_input1[5] / divider_input2[5];
    Pcon1PUfilt = -1.0 * divider_input1[7] / divider_input2[7];
    Qcon1PUfilt = -1.0 * divider_input1[8] / divider_input2[8];
    VacPUfilt = VacRMSfilt[idx] / VLLcon1;
    //VdcBusPU = divider_input1[6] / divider_input2[6];
    
    Kf = 2.0 * (*T_VSC) / *DELT;
    double VdcBusPUfilt;
    LPF_Execute(Kf, *G_VSC, VdcBusPU, &output_his_LPF_VSC[idx], &input_his_LPF_VSC[idx], &VdcBusPUfilt);

    //printf("Thread %d: VdcBusPUfilt[%d] = %f\n", idx, idx, VdcBusPUfilt);

    //PLL and PQ
    double thetaph, thetaLL;
    double VdPUfilt, VqPUfilt, IdPUfilt, IqPUfilt;
    double VconD;
    
    double abc0[2], abc1[2], abc2[2], theta[2];
    double Kf_PLL[4], input_realpole[4], G_PLL[4];
    double input_angleresolver[2];

    PLL3Phase_Execute(&componentPLL3Phase[idx], *Kp_PLL, *Ti_PLL, *Base_Vpll, *fbase, *fmax, *fmin, *DELT, Vcon1ABpu, Vcon1BCpu, Vcon1CApu);
    
    input_angleresolver[0] = componentPLL3Phase[idx].output - PI / 6.0;
    input_angleresolver[1] = componentPLL3Phase[idx].output;

    angle_resolver_Execute(input_angleresolver[0], &thetaph);
    angle_resolver_Execute(input_angleresolver[1], &thetaLL);

    abc0[0] = Vcon1ABpu;
    abc1[0] = Vcon1BCpu;
    abc2[0] = Vcon1CApu;
    abc0[1] = Icon1Apu;
    abc1[1] = Icon1Bpu;
    abc2[1] = Icon1Cpu;
    theta[0] = thetaLL;
    theta[1] = thetaph;

    /*
    if(idx == 0)
    {
        printf("Thread %d: thetaPh[%d] = %f\n", idx, idx, thetaph);
        printf("Thread %d: thetaLL[%d] = %f\n", idx, idx, thetaLL);
    }
    */

    ABC_DQO_Execute(abc0[0], abc1[0], abc2[0], theta[0], &input_realpole[0], &input_realpole[1]);
    ABC_DQO_Execute(abc0[1], abc1[1], abc2[1], theta[1], &input_realpole[2], &input_realpole[3]);

    VconD = input_realpole[0];

    Kf_PLL[0] = Kf_PLL[1] = 2.0 * (*Ti_realpole_L) / (*DELT);
    Kf_PLL[2] = Kf_PLL[3] = 2.0 * (*Ti_realpole_R) / (*DELT);
    G_PLL[0] = G_PLL[1] = G_PLL[2] = G_PLL[3] = *Kp_realpole;

    double output_LPF_PLL[4];
    for (int i = 0; i < 4; i++)
    {
        LPF_Execute(Kf_PLL[i], G_PLL[i], input_realpole[i], &output_his_LPF_PLL[idx * 4 + i], &input_his_LPF_PLL[idx * 4 + i], &output_LPF_PLL[i]);
    }

    VdPUfilt = output_LPF_PLL[0];
    VqPUfilt = -output_LPF_PLL[1];
    IdPUfilt = output_LPF_PLL[2];
    IqPUfilt = -output_LPF_PLL[3];

    /*
    if(idx == 0)
    {
        printf("Thread %d: VdPUfilt[%d] = %f\n", idx, idx, VdPUfilt);
        printf("Thread %d: VqPUfilt[%d] = %f\n", idx, idx, VqPUfilt);
        printf("Thread %d: IdPUfilt[%d] = %f\n", idx, idx, IdPUfilt);
        printf("Thread %d: IqPUfilt[%d] = %f\n", idx, idx, IqPUfilt);
    }
    */
    
    //Voltage Loop
    double signal_switch_input1[2], signal_switch_input2[2];
    int signal_switch_signal[2], ctrl[2];
    double PI_model_in[2], PI_model_RST[2], PI_model_FRZ[2], PI_model_output[2], PI_model_MAX[2], PI_model_MIN[2], PI_model_T[2], PI_model_Kp[2];
    double abs_input[2], abs_output[2];

    double ErrVdcBusPU;
    double ErrVacPU;

    double IdPURef, IqPURef;

    /*
    if(idx == 0)
    {
        printf("Thread %d: VdcBusPUfilt[%d] = %f\n", idx, idx, VdcBusPUfilt);
    }
    */

    signal_switch_input1[0] = *VdcPURef;
    signal_switch_input1[1] = Qcon1PUfilt - *QPURef;
    signal_switch_input2[0] = VdcBusPUfilt;
    signal_switch_input2[1] = *VacPURef - VacPUfilt;

    signal_switch_signal[0] = Dblk1[idx];
    signal_switch_signal[1] = *VacOrQ;
    ctrl[0] = 1;
    ctrl[1] = 0;

    double signal_switch_output[2];
    signal_switch_Execute(signal_switch_input1[0], signal_switch_input2[0], ctrl[0], signal_switch_signal[0], &signal_switch_output[0]);
    signal_switch_Execute(signal_switch_input1[1], signal_switch_input2[1], ctrl[1], signal_switch_signal[1], &signal_switch_output[1]);

    double output_rate_limiter_Voltage;
    rate_limiter_Execute(signal_switch_output[0], *DELT, *ULIM, *LLIM, &output_his_rate_limiter_Voltage[idx], &output_rate_limiter_Voltage);

    ErrVdcBusPU = output_rate_limiter_Voltage - VdcBusPUfilt;
    ErrVacPU = signal_switch_output[1] * StartPI1;

    /*
    if(idx == 0)
    {
        printf("Thread %d: ErrVdcBusPU[%d] = %f\n", idx, idx, ErrVdcBusPU);
        //printf("Thread %d: ErrVacPU[%d] = %f\n", idx, idx, ErrVacPU);
    }
    */

    PI_model_in[0] = ErrVdcBusPU;
    PI_model_in[1] = ErrVacPU;
    PI_model_RST[0] = PI_model_RST[1] = ResetPI1;
    PI_model_FRZ[0] = PI_model_FRZ[1] = *FrzPI;
    PI_model_Kp[0] = *Kp;
    PI_model_Kp[1] = *Kp2;
    PI_model_MAX[0] = *MAXOUT;
    PI_model_MAX[1] = *MAXOUT2;
    PI_model_MIN[0] = *MINOUT;
    PI_model_MIN[1] = *MINOUT2;
    PI_model_T[0] = *PI_T1;
    PI_model_T[1] = *PI_T2;

    for(int i = 0; i < 2; i++)
    {
        PI_model_Execute(&componentPI_model_Voltage[idx * 2 + i], PI_model_in[i], PI_model_RST[i], PI_model_FRZ[i], PI_model_MAX[i], PI_model_MIN[i], PI_model_T[i], PI_model_Kp[i], *DELT);
    }
   
    double output_M_rec_polar_voltage;
    double output_P_rec_polar_voltage;
    rec_polar_Execute(componentPI_model_Voltage[idx * 2 + 0].PI_OUT, componentPI_model_Voltage[idx * 2 + 1].PI_OUT, &output_M_rec_polar_voltage, &output_P_rec_polar_voltage);

    double Limit_MAX;
    if (fabs(*IdMax) >= fabs(*IdMin))
        Limit_MAX = *IdMax;
    else
        Limit_MAX = *IdMin;

    double output_limiter_voltage;
    limiter_Execute(output_M_rec_polar_voltage, Limit_MAX, 0.0, &output_limiter_voltage);

    polar_rec_Execute(output_limiter_voltage, output_P_rec_polar_voltage, &IdPURef, &IqPURef);

    /*
    if(idx == 0)
    {
        printf("Thread %d: IdPURef[%d] = %f\n", idx, idx, IdPURef);
        printf("Thread %d: IqPURef[%d] = %f\n", idx, idx, IqPURef);
    }
    */

    //ErrVdcBusPU_his[idx] = ErrVdcBusPU;
    //ErrVacPU_his[idx] = ErrVacPU;
    //PI_model_output_his1[idx] = componentPI_model_Voltage[idx * 2 + 0].PI_OUT;
    //PI_model_output_his2[idx] = componentPI_model_Voltage[idx * 2 + 1].PI_OUT;

    /*
    if(idx == 0)
    {
        printf("Thread %d: PI_model_output_his1[%d] = %f\n", idx, idx, PI_model_output_his1[0]);
        printf("Thread %d: PI_model_output_his2[%d] = %f\n", idx, idx, PI_model_output_his2[1]);
    }
    */

    //Current Loop
    double thetaCon1;
    double ModWvAx, ModWvBx, ModWvCx;
    
    double Rec_x, Rec_y;
    double limit_input;
    double Mindex;

    double PI_in[2], PI_Max[2], PI_Min[2], PI_Kp[2], PI_Ti[2];

    PI_in[0] = IdPURef - IdPUfilt;
    PI_in[1] = IqPURef - IqPUfilt;
    PI_Max[0] = *MaxLimit1_current;
    PI_Max[1] = *MaxLimit2_current;
    PI_Min[0] = *MinLimit1_current;
    PI_Min[1] = *MinLimit2_current;
    PI_Kp[0] = *Kp1_current;
    PI_Kp[1] = *Kp2_current;
    PI_Ti[0] = *Ti1_current;
    PI_Ti[1] = *Ti2_current;

    for(int i = 0; i < 2; i++)
    {
        PI_model_Execute(&componentPI_model_Current[idx * 2 + i], PI_in[i], ResetPI1, *FrzPI, PI_Max[i], PI_Min[i], PI_Ti[i], PI_Kp[i], *DELT);
    }

    Rec_x = VdPUfilt + IqPUfilt * wLeqPU - componentPI_model_Current[idx * 2 + 0].PI_OUT;
    Rec_y = VqPUfilt - IdPUfilt * wLeqPU - componentPI_model_Current[idx * 2 + 1].PI_OUT;

    double output_M_rec_polar_current;
    double output_P_rec_polar_current;
    rec_polar_Execute(Rec_x, Rec_y, &output_M_rec_polar_current, &output_P_rec_polar_current);

    limit_input = output_M_rec_polar_current * Vphcon1PK / (0.5 * (*VdcRated));
    thetaCon1 = output_P_rec_polar_current + thetaph;

    /*
    if(idx == 0)
    {
        printf("Thread %d: thetaCon1[%d] = %f\n", idx, idx, thetaCon1);
    }
    */

    limiter_Execute(limit_input, *Mmax_current, *Mmin_current, &Mindex);

    double output_X_polar_rec_current, output_Y_polar_rec_current;
    polar_rec_Execute(Mindex, output_P_rec_polar_current, &output_X_polar_rec_current, &output_Y_polar_rec_current);

    DQ0_ABC_Execute(0.0, output_X_polar_rec_current, output_Y_polar_rec_current, thetaph, &ModWvAx, &ModWvBx, &ModWvCx);

    /*
    if(idx == 0)
    {
        printf("Thread %d: ModWvAx[%d] = %f\n", idx, idx, ModWvAx);
        printf("Thread %d: ModWvBx[%d] = %f\n", idx, idx, ModWvBx);
        printf("Thread %d: ModWvCx[%d] = %f\n", idx, idx, ModWvCx);
    }
    */

    //3rd Harmonic
    double ModWvA, ModWvB, ModWvC;
    double signal_switch_input3, signal_switch_input4;

    double output_angle_resolver_threeHarmonic;
    angle_resolver_Execute(3 * thetaCon1, &output_angle_resolver_threeHarmonic);

    double output_signal_switch_threeHarmonic;
    signal_switch_input3 = 0.0;
    signal_switch_input4 = 0.1667 * Mindex * cos(output_angle_resolver_threeHarmonic);
    signal_switch_Execute(signal_switch_input3, signal_switch_input4, *Ctrl, *Har3ctrl, &output_signal_switch_threeHarmonic);

    ModWvA = ModWvAx - output_signal_switch_threeHarmonic;
    ModWvB = ModWvBx - output_signal_switch_threeHarmonic;
    ModWvC = ModWvCx - output_signal_switch_threeHarmonic;

    /*
    if(idx == 0)
    {
        printf("Thread %d: ModWvA[%d] = %f\n", idx, idx, ModWvA);
        printf("Thread %d: ModWvB[%d] = %f\n", idx, idx, ModWvB);
        printf("Thread %d: ModWvC[%d] = %f\n", idx, idx, ModWvC);
    }
    */

    //Boost
    double state;
    double Dbst;
    double input_PI;
    double output_Vmppt;

    double Kf_Boost[3], G_Boost[3], input_real_pole[3];
    Kf_Boost[0] = 2.0 * (*T1_Boost) / (*DELT);
    Kf_Boost[1] = 2.0 * (*T3_Boost) / (*DELT);
    Kf_Boost[2] = 2.0 * (*T2_Boost) / (*DELT);

    G_Boost[0] = *G1_Boost;
    G_Boost[1] = *G3_Boost;
    G_Boost[2] = *G2_Boost;

    input_real_pole[0] = VpvBus[idx];
    input_real_pole[1] = (VpvBus[idx]) / (*VpvRated);
    input_real_pole[2] = Ipv[idx];

    //LPF
    double output_LPF_boost[2];
    double output_LPF2_boost;
    
    for(int i = 0; i < 2; i++)
    {
        LPF_Execute(Kf_Boost[i], G_Boost[i], input_real_pole[i], &output_his_LPF_boost[idx * 2 + i], &input_his_LPF_boost[idx * 2 + i], &output_LPF_boost[i]);
    }
    
    LPF_Execute(Kf_Boost[2], G_Boost[2], input_real_pole[2], &output_his_LPF2_boost[idx], &input_his_LPF2_boost[idx], &output_LPF2_boost);

    //MPPT
    if(*TIME == 0)
    {
        componentMPPT[idx].output_SampNum = int(*Samp * 0.001 / *DELT) + 1;
        componentMPPT[idx].output_StrtupNum = int(*Strtup / *Samp * 1000) + 1;
        componentMPPT[idx].SampleCount = componentMPPT[idx].output_SampNum;
        componentMPPT[idx].StrtupCount = componentMPPT[idx].output_StrtupNum;
    }

    MPPT_Execute(&componentMPPT[idx], Dblk1[idx], output_LPF_boost[0], output_LPF2_boost, *Sd, *Vstep);

    output_Vmppt = componentMPPT[idx].output_Vmppt;

    double output_signal_switch_boost;
    signal_switch_Execute(*signal_A, output_Vmppt, *VdcCon, *signal, &output_signal_switch_boost);

    input_PI = (- output_signal_switch_boost / (*VpvRated) + output_LPF_boost[1]) * StartPI2;
    
    PI_model_Execute(&componentPI_model_boost[idx], input_PI, ResetPI2, *FrzPI, *Max_Boost, *Min_Boost, *Ti_Boost, *Kp_Boost, *DELT);

    Dbst = componentPI_model_boost[idx].PI_OUT;
    state = componentMPPT[idx].output_State;
    state_his[idx] = state;

    /*
    if(idx == 0)
    {
        printf("Thread %d: Dbst[%d] = %f\n", idx, idx, Dbst);
    }
    */

    //PWM
    double H_OFF[7], H_ON[7];
    double L_OFF[7], L_ON[7];
    int Dblk[7];
    double input_freq[2], max_pwm[2], min_pwm[2];

    double ctrl1, ctrl2;
    input_freq[0] = *fsw1;
    input_freq[1] = *fsw2;
    max_pwm[0] = *max1_pwm;
    max_pwm[1] = *max2_pwm;
    min_pwm[0] = *min1_pwm;
    min_pwm[1] = *min2_pwm;
    
    Sawtooth_Generator_Execute(input_freq[0], max_pwm[0], min_pwm[0], *DELT, &Count_i_Sawtooth_Generator[idx * 2 + 0], &ctrl1);
    Sawtooth_Generator_Execute(input_freq[1], max_pwm[1], min_pwm[1], *DELT, &Count_i_Sawtooth_Generator[idx * 2 + 1], &ctrl2);

    /*
    if(idx == 0)
    {
        printf("Thread %d: ctrl1[%d] = %f\n", idx, idx, ctrl1);
        printf("Thread %d: ctrl2[%d] = %f\n", idx, idx, ctrl2);
    }
    */

    H_OFF[0] = L_ON[0] = ctrl1;
    H_OFF[1] = L_ON[1] = ModWvA;
    H_OFF[2] = L_ON[2] = ctrl1;
    H_OFF[3] = L_ON[3] = ModWvB;
    H_OFF[4] = L_ON[4] = ctrl1;
    H_OFF[5] = L_ON[5] = ModWvC;
    H_OFF[6] = L_ON[6] = ctrl2;
    L_OFF[0] = H_ON[0] = ModWvA;
    L_OFF[1] = H_ON[1] = ctrl1;
    L_OFF[2] = H_ON[2] = ModWvB;
    L_OFF[3] = H_ON[3] = ctrl1;
    L_OFF[4] = H_ON[4] = ModWvC;
    L_OFF[5] = H_ON[5] = ctrl1;
    L_OFF[6] = H_ON[6] = Dbst;

    for (int i = 0; i < 6; i++)
    {
        Dblk[i] = Dblk1[idx];
    }
    Dblk[6] = Dblk2;
    
    Interpolated_Firing_Pulses_Execute(Dblk[0], H_OFF[0], L_OFF[0], H_ON[0], L_ON[0], *DELT, &fx_his_Interpolated_Firing_Pulses[idx * 7 + 0], &Ta_1[idx], &Ta_1_timetag[idx]);
    //Ta_1_timetag[idx] = thetaph;
    Interpolated_Firing_Pulses_Execute(Dblk[1], H_OFF[1], L_OFF[1], H_ON[1], L_ON[1], *DELT, &fx_his_Interpolated_Firing_Pulses[idx * 7 + 1], &Ta_4[idx], &Ta_4_timetag[idx]);
    //Ta_4_timetag[idx] = thetaLL;
    Interpolated_Firing_Pulses_Execute(Dblk[2], H_OFF[2], L_OFF[2], H_ON[2], L_ON[2], *DELT, &fx_his_Interpolated_Firing_Pulses[idx * 7 + 2], &Tb_3[idx], &Tb_3_timetag[idx]);
    Interpolated_Firing_Pulses_Execute(Dblk[3], H_OFF[3], L_OFF[3], H_ON[3], L_ON[3], *DELT, &fx_his_Interpolated_Firing_Pulses[idx * 7 + 3], &Tb_6[idx], &Tb_6_timetag[idx]);
    Interpolated_Firing_Pulses_Execute(Dblk[4], H_OFF[4], L_OFF[4], H_ON[4], L_ON[4], *DELT, &fx_his_Interpolated_Firing_Pulses[idx * 7 + 4], &Tc_5[idx], &Tc_5_timetag[idx]);
    Interpolated_Firing_Pulses_Execute(Dblk[5], H_OFF[5], L_OFF[5], H_ON[5], L_ON[5], *DELT, &fx_his_Interpolated_Firing_Pulses[idx * 7 + 5], &Tc_2[idx], &Tc_2_timetag[idx]);
    Interpolated_Firing_Pulses_Execute(Dblk[6], H_OFF[6], L_OFF[6], H_ON[6], L_ON[6], *DELT, &fx_his_Interpolated_Firing_Pulses[idx * 7 + 6], &MPPT_ctrl[idx], &MPPT_ctrl_timetag[idx]);

    /*
    if(idx == 0)
    {
        printf("Thread %d: Ta_1_timetag[%d] = %f\n", idx, idx, Ta_1_timetag[idx]);
        printf("Thread %d: Ta_4_timetag[%d] = %f\n", idx, idx, Ta_4_timetag[idx]);
        printf("Thread %d: Tb_3_timetag[%d] = %f\n", idx, idx, Tb_3_timetag[idx]);
        printf("Thread %d: Tb_6_timetag[%d] = %f\n", idx, idx, Tb_6_timetag[idx]);
        printf("Thread %d: Tc_2_timetag[%d] = %f\n", idx, idx, Tc_2_timetag[idx]);
        printf("Thread %d: Tc_5_timetag[%d] = %f\n", idx, idx, Tc_5_timetag[idx]);
        printf("Thread %d: MPPT_ctrl_timetag[%d] = %f\n", idx, idx, MPPT_ctrl_timetag[idx]);
    }
    */

}

__global__ void SW_AC_DC(int* SW_PV, int* SW_Chopper, int* SW_Char, int* MPPT_ctrl, int* Ta_1, int* Tb_3, int* Tc_5, int* Ta_4, int* Tb_6, int* Tc_2, int* SW_PVunit, double* VdcPreCh,
    int* SW_DC, int* SW_AC, double* Vchar)
{
    int idx = threadIdx.x;
            
    SW_DC[idx * 4 + 0] = SW_PV[idx];
    SW_DC[idx * 4 + 1] = SW_Chopper[idx];
    SW_DC[idx * 4 + 2] = SW_Char[idx];
    SW_DC[idx * 4 + 3] = MPPT_ctrl[idx];

    SW_AC[idx * 7 + 0] = Ta_1[idx];
    SW_AC[idx * 7 + 1] = Tb_3[idx];
    SW_AC[idx * 7 + 2] = Tc_5[idx];
    SW_AC[idx * 7 + 3] = Ta_4[idx];
    SW_AC[idx * 7 + 4] = Tb_6[idx];
    SW_AC[idx * 7 + 5] = Tc_2[idx]; 
    SW_AC[idx * 7 + 6] = SW_PVunit[idx];

    Vchar[idx] = VdcPreCh[idx];

}

__global__ void pvmodel_pre_dc(double* _DELT, double* _PVarray_pars, double* _Boost_pars, double* _CPB_pars, double* _VSC_pars,
    Switch* Bre_PV, Switch* Bre_Chop, Switch* Bre_Char, Switch* Boost_IGBT, Switch* Boost_Diode, Capacitor* C_PV, Capacitor* C_Boost, Capacitor* C_VSC, Inductor* L_Boost,
    double* Npar, double* Nser, double* Isc1_PV, double* Voc1_PV, double* par_PV_c1, double* par_PV_b1, double* V_PV_V_his, double* I_PV_A, double* I_PV_kA, double* STATUS_his,
    double** d_Yg2_d, double** d_J_branch2_d)
{
	int idx = threadIdx.x;
	
    //Initializing parameter
	//PV
	//PV Breaker [intial status: OPEN]
	double Roff_PV = _PVarray_pars[14];
	double Ron_PV = _PVarray_pars[15];
	//Compared with RTDS?
	Bre_PV[idx].gon = 1.0 / Ron_PV;
	Bre_PV[idx].goff = 1.0 / Roff_PV;
	Bre_PV[idx].g = Bre_PV[idx].goff;

	//PV Capacitor
	double CPV = _PVarray_pars[16];

	C_PV[idx].cap_g = 2.0 * CPV * 1.0e-6 / *_DELT;
	C_PV[idx].cap_jh = 0.0;
	C_PV[idx].cap_ib = 0.0;
	C_PV[idx].cap_v = 0.0;

	//PVarray
	double Voc_PV, Isc_PV, Vm_PV, Im_PV;
	double Tref, Sref;
	double INSO_ini, TEMP_ini, Vpv_ini;
	double par_PV_a, par_PV_b, par_PV_c;

	Voc_PV = _PVarray_pars[0];
	Isc_PV = _PVarray_pars[1];
	Vm_PV = _PVarray_pars[2];
	Im_PV = _PVarray_pars[3];
	Npar[idx] = _PVarray_pars[4];
	Nser[idx] = _PVarray_pars[5];
	Tref = _PVarray_pars[6];
	Sref = _PVarray_pars[7];
	INSO_ini = _PVarray_pars[8];
	TEMP_ini = _PVarray_pars[9]; 
	Vpv_ini = _PVarray_pars[10];
	par_PV_a = _PVarray_pars[11];
	par_PV_b = _PVarray_pars[12];
	par_PV_c = _PVarray_pars[13];

	double Im1_PV, Vm1_PV;
	double dTemp, SdSref, dS;

	dTemp = TEMP_ini - Tref;
	SdSref = INSO_ini / Sref;
	dS = SdSref - 1.0;

	//Modified Expression I'sc V'oc I'm V'm
	Isc1_PV[idx] = Isc_PV * SdSref * (par_PV_a * dTemp + 1.0);
	Im1_PV = Im_PV * SdSref * (par_PV_a * dTemp + 1.0);
	Voc1_PV[idx] = Voc_PV * (1.0 - dTemp * par_PV_c) * log(M_E + par_PV_b * dS);
	Vm1_PV = Vm_PV * (1.0 - dTemp * par_PV_c) * log(M_E + par_PV_b * dS);

	//Calculate C1(par_PV_b1) and C2(par_PV_c1)
	par_PV_c1[idx] = (Vm1_PV / Voc1_PV[idx] - 1.0) / log(1.0 - Im1_PV / Isc1_PV[idx]);
	par_PV_b1[idx] = (1 - Im1_PV / Isc1_PV[idx]) * exp(Vm1_PV / (-1.0 * par_PV_c1[idx] * Voc1_PV[idx]));

	V_PV_V_his[idx] = Vpv_ini * 1000.0 / Nser[idx];
	// I = Isc * (1 - C1 * (exp(V / C2 / Voc)-1))
	I_PV_A[idx] = Isc1_PV[idx] * (1.0 - par_PV_b1[idx] * (exp(V_PV_V_his[idx] / par_PV_c1[idx] / Voc1_PV[idx]) - 1.0));

	//Boost
	//Boost Inductor [unit: H]
	double LBoost = _Boost_pars[0];

	L_Boost[idx].ind_g = *_DELT / (2.0 * LBoost);
	L_Boost[idx].ind_jh = 0.0;
	L_Boost[idx].ind_ib = 0.0;
	L_Boost[idx].ind_v = 0.0;

	//Boost Capacitor [unit: uF]
	double CBoost = _Boost_pars[1];

	C_Boost[idx].cap_g = 2.0 * CBoost * 1.0e-6 / *_DELT;
	C_Boost[idx].cap_jh = 0.0;
	C_Boost[idx].cap_ib = 0.0;
	C_Boost[idx].cap_v = 0.0;

	//IGBT [intial status: OFF]
	double Ron_IGBT = _Boost_pars[2];
	double Roff_IGBT = _Boost_pars[3];

	Boost_IGBT[idx].gon = 1.0 / Ron_IGBT;
	Boost_IGBT[idx].goff = 1.0 / Roff_IGBT;
	Boost_IGBT[idx].g = Boost_IGBT[idx].goff;

	//Diode [intial status: OFF]
	double Ron_Diode = _Boost_pars[4];
	double Roff_Diode = _Boost_pars[5];

	Boost_Diode[idx].gon = 1.0 / Ron_Diode;
	Boost_Diode[idx].goff = 1.0 / Roff_Diode;
	Boost_Diode[idx].g = Boost_Diode[idx].goff;

	//Boost STATUS  0:Diode  1:IGBT
	STATUS_his[idx * 2] = 0;
	STATUS_his[idx * 2 + 1] = 0;

	//Chopper Branch [intial status: OPEN]
	double Roff_Chop = _CPB_pars[0];
	double Ron_Chop = _CPB_pars[1];

	Bre_Chop[idx].gon = 1.0 / Ron_Chop;
	Bre_Chop[idx].goff = 1.0 / Roff_Chop;
	Bre_Chop[idx].g =  Bre_Chop[idx].goff;

	//Precharge Branch[intial status: OPEN]
	double Roff_Char = _CPB_pars[2];
	double Ron_Char = _CPB_pars[3];

	Bre_Char[idx].gon = 1.0 / Ron_Char;
	Bre_Char[idx].goff = 1.0 / Roff_Char;
	Bre_Char[idx].g =  Bre_Char[idx].gon;
	
	//ATTETION
	double j_Char = 0.0;

	//VSC
	//DC Bus Capacitor
	double CBus = _VSC_pars[3];

	C_VSC[idx].cap_g = 2.0 * CBus * 1.0e-6 / *_DELT;
	C_VSC[idx].cap_jh = 0.0;
	C_VSC[idx].cap_ib = 0.0;
	C_VSC[idx].cap_v = 0.0;

	//Switch
	double Ron_VSC = _VSC_pars[0];
	double Roff_VSC = _VSC_pars[1];
	double Rsnubber = _VSC_pars[2];

	//Generate Yg
	d_Yg2_d[idx][0 + 0 * BRANCH_NUM_DC] = C_VSC[idx].cap_g; //[NO CHANGE]

	//Precharge and Chopper
	d_Yg2_d[idx][1 + 1 * BRANCH_NUM_DC] = Bre_Char[idx].g;
	d_Yg2_d[idx][2 + 2 * BRANCH_NUM_DC] = Bre_Chop[idx].g;

	//Boost
	d_Yg2_d[idx][3 + 3 * BRANCH_NUM_DC] = C_Boost[idx].cap_g;
	d_Yg2_d[idx][4 + 4 * BRANCH_NUM_DC] = Boost_Diode[idx].g;
	d_Yg2_d[idx][5 + 5 * BRANCH_NUM_DC] = Boost_IGBT[idx].g;
	d_Yg2_d[idx][6 + 6 * BRANCH_NUM_DC] = L_Boost[idx].ind_g;

	//PV array
	d_Yg2_d[idx][7 + 7 * BRANCH_NUM_DC] = Bre_PV[idx].g;
	d_Yg2_d[idx][8 + 8 * BRANCH_NUM_DC] = C_PV[idx].cap_g;

	
    //Generate J_branch2
	d_J_branch2_d[idx][0] = C_VSC[idx].cap_jh;
	d_J_branch2_d[idx][1] = j_Char; //Depend on Signal
    d_J_branch2_d[idx][2] = 0;
	d_J_branch2_d[idx][3] = C_Boost[idx].cap_jh;
    d_J_branch2_d[idx][4] = 0;
    d_J_branch2_d[idx][5] = 0;
	d_J_branch2_d[idx][6] = L_Boost[idx].ind_jh;
    d_J_branch2_d[idx][7] = 0;
	//ATTENTION //Please Check
	I_PV_kA[idx] = I_PV_A[idx] * Npar[idx] * 0.001;
	d_J_branch2_d[idx][8] = C_PV[idx].cap_jh + I_PV_kA[idx];

}

__global__ void pvmodel_pre_ac(double* _DELT, double* _VSC_pars, double* _LC_pars, double* _PVunit_Bre_pars, double* _Trans_pars, 
    Switch* Bre_PVunit, VSC* VSC_average, Capacitor* C_filt, Inductor* L_filt,
    double* G_filt_R1, double* G_filt_R2, double* d_G_trans, double* d_parB_trans, double* d_J_trans_A_PVunit_d, double* d_J_trans_B_PVunit_d, double* d_J_trans_C_PVunit_d,
    double** d_Yg1_d, double** d_J_branch1_d)
{
	int idx = threadIdx.x; 

	//Initializing parameter
	//VSC
	//Switch
	double Ron_VSC = _VSC_pars[0];
	double Roff_VSC = _VSC_pars[1];
	double Rsnubber = _VSC_pars[2];
	
	VSC_average[idx].gon_vsc = 1.0 / Ron_VSC;
	VSC_average[idx].goff_vsc = 1.0 / Roff_VSC;
	VSC_average[idx].g_snubber = 1.0 / Rsnubber;
	VSC_average[idx].g_vsc = VSC_average[idx].g_snubber;

	//LC filter
	double L_LC = _LC_pars[0];
	double R1_LC = _LC_pars[1];
	double C_LC = _LC_pars[2];
	double R2_LC = _LC_pars[3];

	G_filt_R1[idx] = 1.0 / R1_LC;
	G_filt_R2[idx] = 1.0 / R2_LC;

	//Inductor
	for (int i = 0; i < 3; i++)
	{
		L_filt[idx * 3 + i].ind_g = *_DELT / (2.0 * L_LC);
		L_filt[idx * 3 + i].ind_v = 0.0;
		L_filt[idx * 3 + i].ind_ib = 0.0;
		L_filt[idx * 3 + i].ind_jh = 0.0;
	}
	
	//Capacitor
	for (int i = 0; i < 3; i++)
	{
		C_filt[idx * 3 + i].cap_g = 2.0 * C_LC * 1.0e-6 / *_DELT;
		C_filt[idx * 3 + i].cap_v = 0.0;
		C_filt[idx * 3 + i].cap_ib = 0.0;
		C_filt[idx * 3 + i].cap_jh = 0.0;
	}

	//PV unit breaker
	double Roff_PVunit = _PVunit_Bre_pars[0];
	double Ron_PVunit = _PVunit_Bre_pars[1];
	
	Bre_PVunit[idx].gon = 1.0 / Ron_PVunit;
	Bre_PVunit[idx].goff = 1.0 / Roff_PVunit;
	Bre_PVunit[idx].g = Bre_PVunit[idx].goff;

	//Transformer
	double Sb_trans = _Trans_pars[0];
	double fb_trans = _Trans_pars[1];
	double XL_trans = _Trans_pars[2];
	double Z_trans = _Trans_pars[3];
	double Im_trans = _Trans_pars[4];
	double V1_trans = _Trans_pars[5];
	double V2_trans = _Trans_pars[6];

	double N_trans = V1_trans / V2_trans;
	double Zb1_trans = V1_trans * V1_trans / Sb_trans;
	double Zb2_trans = V2_trans * V2_trans / Sb_trans;

	double L1_trans = (XL_trans / 2) * Zb1_trans / (2 * PI * fb_trans);
	double L2_trans = L1_trans / N_trans / N_trans;
	double Lm_trans = 1.0 / Im_trans * Zb1_trans / (2 * PI * fb_trans); 

	double R1_trans = Z_trans / 2 * Zb1_trans;
	double R2_trans = Z_trans / 2 * Zb2_trans;

	double R_trans[2][2];
	double L_trans[2][2];

	R_trans[0][0] = R1_trans;
    R_trans[0][1] = 0;
    R_trans[1][0] = 0;
	R_trans[1][1] = R2_trans;

	L_trans[0][0] = L1_trans + Lm_trans;
	L_trans[0][1] = Lm_trans / N_trans;
	L_trans[1][0] = Lm_trans / N_trans;
	L_trans[1][1] = L2_trans + Lm_trans / N_trans / N_trans;

	double parA_trans[2][2];
    double parB_trans[2][2];
    double G_trans[2][2];

	for(int i = 0; i < 2; i++)
    {
        for(int j = 0; j < 2; j++)
        {
            parA_trans[i][j] = R_trans[i][j] + 2 / (*_DELT) * L_trans[i][j];
            parB_trans[i][j] = -R_trans[i][j] +  2 / (*_DELT) * L_trans[i][j];
        }
    }

    double inva = parA_trans[0][0];
    double invb = parA_trans[0][1];
    double invc = parA_trans[1][0];
    double invd = parA_trans[1][1];
    double invdet = inva * invd - invb * invc;
    double inva1 = invd / invdet;
    double invb1 = -invb / invdet;
    double invc1 = -invc / invdet;
    double invd1 = inva / invdet;

    G_trans[0][0] = inva1;
    G_trans[0][1] = invb1;
    G_trans[1][0] = invc1;
    G_trans[1][1] = invd1;

	for (int i = 0; i < 2; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            d_G_trans[idx * 4 + i + j * 2] = G_trans[i][j];
        }
    }

    for (int i = 0; i < 2; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            d_parB_trans[idx * 4 + i + j * 2] = parB_trans[i][j];
        }
    }
	
    for (int i = 0; i < 2; i++)
    {
        d_J_trans_A_PVunit_d[idx * 2 + i] = 0;
        d_J_trans_B_PVunit_d[idx * 2 + i] = 0;
        d_J_trans_C_PVunit_d[idx * 2 + i] = 0;
    }

	//Generate Yg
	//Transformer
	for(int i = 0; i < 3; i++)
	{
		d_Yg1_d[idx][i * 6 + i * 6 * BRANCH_NUM_AC] = d_G_trans[idx * 4 + 1 + 1 * 2];
		d_Yg1_d[idx][i * 6 + 1 + (i * 6 + 1) * BRANCH_NUM_AC] = -d_G_trans[idx * 4 + 0 + 1 * 2];
		d_Yg1_d[idx][i * 6 + 2 + (i * 6 + 2) * BRANCH_NUM_AC] = d_G_trans[idx * 4 + 0 + 1 * 2];
		d_Yg1_d[idx][i * 6 + 3 + (i * 6 + 3) * BRANCH_NUM_AC] = d_G_trans[idx * 4 + 0 + 1 * 2];
		d_Yg1_d[idx][i * 6 + 4 + (i * 6 + 4) * BRANCH_NUM_AC] = -d_G_trans[idx * 4 + 0 + 1 * 2];
		d_Yg1_d[idx][i * 6 + 5 + (i * 6 + 5) * BRANCH_NUM_AC] = d_G_trans[idx * 4 + 0 + 0 * 2];
	}
	
	//PVunit Breaker
	d_Yg1_d[idx][18 + 18 * BRANCH_NUM_AC] = Bre_PVunit[idx].g;
	d_Yg1_d[idx][19 + 19 * BRANCH_NUM_AC] = Bre_PVunit[idx].g;
	d_Yg1_d[idx][20 + 20 * BRANCH_NUM_AC] = Bre_PVunit[idx].g;
	
	//LC filter [NO CHANGE]
	for(int i = 0; i < 3; i++)
	{
		d_Yg1_d[idx][i + 21 + (i + 21) * BRANCH_NUM_AC] = G_filt_R2[idx];
		d_Yg1_d[idx][i + 24 + (i + 24) * BRANCH_NUM_AC] = C_filt[idx * 3 + i].cap_g;
		d_Yg1_d[idx][i + 27 + (i + 27) * BRANCH_NUM_AC] = G_filt_R1[idx];
		d_Yg1_d[idx][i + 30 + (i + 30) * BRANCH_NUM_AC] = L_filt[idx * 3 + i].ind_g;
	}

	//VSC [initial status: off] [CHANGE]
	for(int i = 33; i < 36; i++)
	{
		d_Yg1_d[idx][i + i * BRANCH_NUM_AC] = VSC_average[idx].g_vsc;
	}

	//Generate J_branch
    for(int i = 0; i < BRANCH_NUM_AC; i++)
    {
        d_J_branch1_d[idx][i] = 0;
    }
	
	d_J_branch1_d[idx][0] = d_J_trans_A_PVunit_d[idx * 2 + 1];
	d_J_branch1_d[idx][5] = d_J_trans_A_PVunit_d[idx * 2 + 0];
	d_J_branch1_d[idx][6] = d_J_trans_B_PVunit_d[idx * 2 + 1];
	d_J_branch1_d[idx][11] = d_J_trans_B_PVunit_d[idx * 2 + 0];
	d_J_branch1_d[idx][12] = d_J_trans_C_PVunit_d[idx * 2 + 1];
	d_J_branch1_d[idx][17] = d_J_trans_C_PVunit_d[idx * 2 + 0];

	d_J_branch1_d[idx][24] = C_filt[idx * 3 + 0].cap_jh;
	d_J_branch1_d[idx][25] = C_filt[idx * 3 + 1].cap_jh;
	d_J_branch1_d[idx][26] = C_filt[idx * 3 + 2].cap_jh;

	d_J_branch1_d[idx][30] = L_filt[idx * 3 + 0].ind_jh;
	d_J_branch1_d[idx][31] = L_filt[idx * 3 + 1].ind_jh;
	d_J_branch1_d[idx][32] = L_filt[idx * 3 + 2].ind_jh;

}

__global__ void MatrixPartitioning_Y11(double **d_Y1_d, double **d_Y11_d) {
    int matrixIdx = blockIdx.x; 
    int idx = threadIdx.x;

    int i = idx % 3; // 计算当前线程负责的行
    int j = idx / 3; // 计算当前线程负责的列

    if (i < 3 && j < 3) {
        d_Y11_d[matrixIdx][i + j * 3] = d_Y1_d[matrixIdx][i + j * NODE_NUM_AC];
    }
}

__global__ void MatrixPartitioning_Y12(double **d_Y1_d, double **d_Y12_d) {
    int matrixIdx = blockIdx.x; 
    int idx = threadIdx.x;

    int i = idx % 3;
    int j = idx / 3;

    if (i < 3 && j < (NODE_NUM_AC - 3)) {
        d_Y12_d[matrixIdx][i + j * 3] = d_Y1_d[matrixIdx][i + (j + 3) * NODE_NUM_AC];
    }
}

__global__ void MatrixPartitioning_Y21(double **d_Y1_d, double **d_Y21_d) {
    int matrixIdx = blockIdx.x; 
    int idx = threadIdx.x;

    int i = idx % (NODE_NUM_AC - 3);
    int j = idx / (NODE_NUM_AC - 3);

    if (i < (NODE_NUM_AC - 3) && j < 3) {
        d_Y21_d[matrixIdx][i + j * (NODE_NUM_AC - 3)] = d_Y1_d[matrixIdx][i + 3 + j * NODE_NUM_AC];
    }
}

__global__ void MatrixPartitioning_Y22(double **d_Y1_d, double **d_Y22_d) {
    int matrixIdx = blockIdx.x; 
    int idx = threadIdx.x;

    int i = idx % (NODE_NUM_AC - 3);
    int j = idx / (NODE_NUM_AC - 3);

    if (i < (NODE_NUM_AC - 3) && j < (NODE_NUM_AC - 3)) {
        d_Y22_d[matrixIdx][i + j * (NODE_NUM_AC - 3)] = d_Y1_d[matrixIdx][i + 3 + (j + 3) * NODE_NUM_AC];
    }
}

__global__ void MatrixPartitioning_Jex(double **d_J1_d, double **d_Jex_d) {
    int matrixIdx = blockIdx.x; 
    int idx = threadIdx.x;

    if (idx < 3) {
        d_Jex_d[matrixIdx][idx] = d_J1_d[matrixIdx][idx];
    }
}

__global__ void MatrixPartitioning_Jin(double **d_J1_d, double **d_Jin_d) {
    int matrixIdx = blockIdx.x; 
    int idx = threadIdx.x;

    if (idx < (NODE_NUM_AC - 3)) {
        d_Jin_d[matrixIdx][idx] = d_J1_d[matrixIdx][idx + 3];
    }
}

__global__ void JYV_Calculation(double **d_Jin_d, double **d_YV_d, double **d_JYV_d)
{
    int idx = threadIdx.x;

    for(int i = 0; i < NODE_NUM_AC - 3; i++)
    {
        d_JYV_d[idx][i] = d_Jin_d[idx][i] - d_YV_d[idx][i];
    }
    
}

__global__ void Yeq_Jeq_Calculation(double **d_Y11_d, double **d_YYY_d, double **d_Yeq_d, double **d_Jex_d, double **d_YYJ_d, double **d_Jeq_d, double *G_PVunit, double *INJ_ABC)
{
    int idx = threadIdx.x;

    for (int i = 0; i < 3; i++) 
    {
        for (int j = 0; j < 3; j++) 
        {
            d_Yeq_d[idx][i + j * 3] = d_Y11_d[idx][i + j * 3] - d_YYY_d[idx][i + j * 3];
        }
    }
        
    for (int i = 0; i < 3; i++) 
    {
        d_Jeq_d[idx][i] = d_Jex_d[idx][i] - d_YYJ_d[idx][i];
    }

    INJ_ABC[idx * 3 + 0] = -d_Jeq_d[idx][0];
	INJ_ABC[idx * 3 + 1] = -d_Jeq_d[idx][1];
    INJ_ABC[idx * 3 + 2] = -d_Jeq_d[idx][2];

	G_PVunit[idx * 6 + 0] = d_Yeq_d[idx][0] + d_Yeq_d[idx][1] + d_Yeq_d[idx][2];
	G_PVunit[idx * 6 + 1] = d_Yeq_d[idx][3] + d_Yeq_d[idx][4] + d_Yeq_d[idx][5];
	G_PVunit[idx * 6 + 2] = d_Yeq_d[idx][6] + d_Yeq_d[idx][7] + d_Yeq_d[idx][8];
	G_PVunit[idx * 6 + 3] = -d_Yeq_d[idx][1];
	G_PVunit[idx * 6 + 4] = -d_Yeq_d[idx][5];
	G_PVunit[idx * 6 + 5] = -d_Yeq_d[idx][2];

}

__global__ void pvmodel_dc(int* _SW_DC, double* _Vchar, double* _idc_source, double** d_V_DC_d,
    Switch* Bre_PV, Switch* Bre_Chop, Switch* Bre_Char, Switch* Boost_IGBT, Switch* Boost_Diode, Capacitor* C_PV, Capacitor* C_Boost, Capacitor* C_VSC, Inductor* L_Boost,
    double* Npar, double* Nser, double* Isc1_PV, double* Voc1_PV, double* par_PV_c1, double* par_PV_b1, double* V_PV_V_his, double* V_PV_V, double* V_PV_kV, double* I_PV_A, double* I_PV_kA, double* STATUS_his,
    double** d_Yg2_d, double** d_J_branch2_d)
{
    int idx = threadIdx.x;
    
    //SIGNAL
	int SW_BRE_PV, SW_BRE_Chop, SW_BRE_Char;
	int SW_Boost;

	SW_BRE_PV = _SW_DC[idx * 4 + 0];
	SW_BRE_Chop = _SW_DC[idx * 4 + 1];
	SW_BRE_Char = _SW_DC[idx * 4 + 2];
	SW_Boost = _SW_DC[idx * 4 + 3];

	//VSC Capacitor
	C_VSC[idx].cap_v = d_V_DC_d[idx][0] - 0;

	//Boost
	C_Boost[idx].cap_v = d_V_DC_d[idx][0] - 0;
	L_Boost[idx].ind_v = d_V_DC_d[idx][2] - d_V_DC_d[idx][1];

	//PV array Capacitor
	C_PV[idx].cap_v = 0 - d_V_DC_d[idx][3];

	//Calculate Branch Current
	//Capacitor
	C_VSC[idx].cap_ib = C_VSC[idx].cap_g * C_VSC[idx].cap_v + C_VSC[idx].cap_jh;
	C_Boost[idx].cap_ib = C_Boost[idx].cap_g * C_Boost[idx].cap_v + C_Boost[idx].cap_jh;
	C_PV[idx].cap_ib = C_PV[idx].cap_g * C_PV[idx].cap_v + C_PV[idx].cap_jh;

	//Inductor
	L_Boost[idx].ind_ib = L_Boost[idx].ind_g * L_Boost[idx].ind_v + L_Boost[idx].ind_jh;

	//Update branch G
	//BREAKER
	Bre_PV[idx].g = SW_BRE_PV * Bre_PV[idx].goff + (1 - SW_BRE_PV) * Bre_PV[idx].gon;
	Bre_Chop[idx].g = SW_BRE_Chop * Bre_Chop[idx].goff + (1 - SW_BRE_Chop) * Bre_Chop[idx].gon;
	Bre_Char[idx].g = SW_BRE_Char * Bre_Char[idx].goff + (1 - SW_BRE_Char) * Bre_Char[idx].gon;

	//Boost
	double iLcheck = 2 * L_Boost[idx].ind_g * L_Boost[idx].ind_v + L_Boost[idx].ind_ib;
	int STATUS[2];

	//IGBT
	if(SW_Boost == 1)
	{
		STATUS[1] = 1;
	}
	else
	{
		STATUS[1] = 0;
	}

	//Diode
	STATUS[0] = 0;
	if(STATUS[1] == 0)
	{
		if((STATUS_his[idx * 2 + 1] == 1) ||(STATUS_his[idx * 2 + 0] == 1 && iLcheck > 0))
		{
			STATUS[0] = 1;
		}
	}

	Boost_IGBT[idx].g = STATUS[1] * Boost_IGBT[idx].gon + (1 - STATUS[1]) * Boost_IGBT[idx].goff;
	Boost_Diode[idx].g = STATUS[0] * Boost_Diode[idx].gon + (1 - STATUS[0]) * Boost_Diode[idx].goff;

	STATUS_his[idx * 2 + 0] = STATUS[0];
	STATUS_his[idx * 2 + 1] = STATUS[1];

	//Gernate Yg
	//VSC
	double idc_source = _idc_source[idx];
	
	//Precharge and Chopper //CHeck Bre_Char
	d_Yg2_d[idx][1 + 1 * BRANCH_NUM_DC] = Bre_Char[idx].g;
	d_Yg2_d[idx][2 + 2 * BRANCH_NUM_DC] = Bre_Chop[idx].g;

	//Boost
	d_Yg2_d[idx][4 + 4 * BRANCH_NUM_DC] = Boost_Diode[idx].g;
	d_Yg2_d[idx][5 + 5 * BRANCH_NUM_DC] = Boost_IGBT[idx].g;
	
	//PV array
	d_Yg2_d[idx][7 + 7 * BRANCH_NUM_DC] = Bre_PV[idx].g;

	//Update History Current
	//Capacitor
	C_VSC[idx].cap_jh = - C_VSC[idx].cap_g * C_VSC[idx].cap_v - C_VSC[idx].cap_ib;
	C_Boost[idx].cap_jh = - C_Boost[idx].cap_g * C_Boost[idx].cap_v - C_Boost[idx].cap_ib;
	C_PV[idx].cap_jh = - C_PV[idx].cap_g * C_PV[idx].cap_v - C_PV[idx].cap_ib;

	//Inductor
	L_Boost[idx].ind_jh = L_Boost[idx].ind_g * L_Boost[idx].ind_v + L_Boost[idx].ind_ib;

	//PV array
	V_PV_kV[idx] = d_V_DC_d[idx][3] - 0;
	V_PV_V[idx] = V_PV_kV[idx] * 1000.0 / Nser[idx];
	V_PV_V_his[idx] = V_PV_V[idx];
	I_PV_A[idx] = Isc1_PV[idx] * (1.0 - par_PV_b1[idx] * (exp(V_PV_V_his[idx] / par_PV_c1[idx] / Voc1_PV[idx]) - 1.0));
	I_PV_kA[idx] =  I_PV_A[idx] * Npar[idx] * 0.001;

	//Precharge Branch
	double j_Char = - _Vchar[idx] * Bre_Char[idx].g;

	//Generate J
	d_J_branch2_d[idx][0] = C_VSC[idx].cap_jh + idc_source;

	d_J_branch2_d[idx][1] = j_Char; //Depend on Signal

	d_J_branch2_d[idx][3] = C_Boost[idx].cap_jh;
	d_J_branch2_d[idx][6] = L_Boost[idx].ind_jh;

	//ATTENTION //Please Check
	d_J_branch2_d[idx][8] = C_PV[idx].cap_jh + I_PV_kA[idx];

}

__global__ void pvmodel_ac(int* _SW_AC, int* _Dblk1, double* _VdcBus, double* _idc_source,
    double** d_Vex_d, double** d_Vin_d, double* I_port_A_PVunit, double* I_port_B_PVunit, double* I_port_C_PVunit,
    Switch* Bre_PVunit, VSC* VSC_average, Capacitor* C_filt, Inductor* L_filt,
    double* G_filt_R1, double* G_filt_R2, double* d_G_trans, double* d_parB_trans, double* d_J_trans_A_PVunit_d, double* d_J_trans_B_PVunit_d, double* d_J_trans_C_PVunit_d,
    double** d_Yg1_d, double** d_J_branch1_d)
{
    int idx = threadIdx.x;
    
    //SIGNAL
	int SW_VSC[6];

	for(int i = 0; i < 6; i++)
	{
		SW_VSC[i] = _SW_AC[idx * 7 + i];
	}

	int SW_PVunit;

	SW_PVunit = _SW_AC[idx * 7 + 6]; 
	
	//Calculate branch voltage and current
	//Transformer
	double V_port_A[2], V_port_B[2], V_port_C[2];

	V_port_A[0] = d_Vin_d[idx][0] - 0;
	V_port_B[0] = d_Vin_d[idx][1] - 0;
	V_port_C[0] = d_Vin_d[idx][2] - 0;

	V_port_A[1] = d_Vex_d[idx][0] - d_Vex_d[idx][1];
	V_port_B[1] = d_Vex_d[idx][1] - d_Vex_d[idx][2];
	V_port_C[1] = d_Vex_d[idx][2] - d_Vex_d[idx][0];

	//LC filter
	C_filt[idx * 3 + 0].cap_v = d_Vin_d[idx][6] - 0;
	C_filt[idx * 3 + 1].cap_v = d_Vin_d[idx][7] - 0;
	C_filt[idx * 3 + 2].cap_v = d_Vin_d[idx][8] - 0;

	L_filt[idx * 3 + 0].ind_v = d_Vin_d[idx][9] - d_Vin_d[idx][12];
	L_filt[idx * 3 + 1].ind_v = d_Vin_d[idx][10] - d_Vin_d[idx][13];
	L_filt[idx * 3 + 2].ind_v = d_Vin_d[idx][11] - d_Vin_d[idx][14];

	//Calculate Branch Current
	//Transformer
	I_port_A_PVunit[idx * 2 + 0] = d_G_trans[idx * 4 + 0] * V_port_A[0] + d_G_trans[idx * 4 + 2] * V_port_A[1] + d_J_trans_A_PVunit_d[idx * 2 + 0];
    I_port_A_PVunit[idx * 2 + 1] = d_G_trans[idx * 4 + 1] * V_port_A[0] + d_G_trans[idx * 4 + 3] * V_port_A[1] + d_J_trans_A_PVunit_d[idx * 2 + 1];
	I_port_B_PVunit[idx * 2 + 0] = d_G_trans[idx * 4 + 0] * V_port_B[0] + d_G_trans[idx * 4 + 2] * V_port_B[1] + d_J_trans_B_PVunit_d[idx * 2 + 0];
    I_port_B_PVunit[idx * 2 + 1] = d_G_trans[idx * 4 + 1] * V_port_B[0] + d_G_trans[idx * 4 + 3] * V_port_B[1] + d_J_trans_B_PVunit_d[idx * 2 + 1];
    I_port_C_PVunit[idx * 2 + 0] = d_G_trans[idx * 4 + 0] * V_port_C[0] + d_G_trans[idx * 4 + 2] * V_port_C[1] + d_J_trans_C_PVunit_d[idx * 2 + 0];
    I_port_C_PVunit[idx * 2 + 1] = d_G_trans[idx * 4 + 1] * V_port_C[0] + d_G_trans[idx * 4 + 3] * V_port_C[1] + d_J_trans_C_PVunit_d[idx * 2 + 1];

	//Capacitor
	for(int i = 0; i < 3; i++)
	{
		C_filt[idx * 3 + i].cap_ib = C_filt[idx * 3 + i].cap_g * C_filt[idx * 3 + i].cap_v + C_filt[idx * 3 + i].cap_jh;
	}

	//Inductor
	for(int i = 0; i < 3; i++)
	{
		L_filt[idx * 3 + i].ind_ib = L_filt[idx * 3 + i].ind_g * L_filt[idx * 3 + i].ind_v + L_filt[idx * 3 + i].ind_jh;
	}

	//Update History Current
	//Transformer
	//d_J_trans_A_array_d
    double V_port_temp_A[2];
    V_port_temp_A[0] = V_port_A[0] + d_parB_trans[idx * 4 + 0] * I_port_A_PVunit[idx * 2 + 0] + d_parB_trans[idx * 4 + 2] * I_port_A_PVunit[idx * 2 + 1];
    V_port_temp_A[1] = V_port_A[1] + d_parB_trans[idx * 4 + 1] * I_port_A_PVunit[idx * 2 + 0] + d_parB_trans[idx * 4 + 3] * I_port_A_PVunit[idx * 2 + 1];
    d_J_trans_A_PVunit_d[idx * 2 + 0] = d_G_trans[idx * 4 + 0] * V_port_temp_A[0] + d_G_trans[idx * 4 + 2] * V_port_temp_A[1];
    d_J_trans_A_PVunit_d[idx * 2 + 1] = d_G_trans[idx * 4 + 1] * V_port_temp_A[0] + d_G_trans[idx * 4 + 3] * V_port_temp_A[1];

    //d_J_trans_B_array_d
    double V_port_temp_B[2];
    V_port_temp_B[0] = V_port_B[0] + d_parB_trans[idx * 4 + 0] * I_port_B_PVunit[idx * 2 + 0] + d_parB_trans[idx * 4 + 2] * I_port_B_PVunit[idx * 2 + 1];
    V_port_temp_B[1] = V_port_B[1] + d_parB_trans[idx * 4 + 1] * I_port_B_PVunit[idx * 2 + 0] + d_parB_trans[idx * 4 + 3] * I_port_B_PVunit[idx * 2 + 1];
    d_J_trans_B_PVunit_d[idx * 2 + 0] = d_G_trans[idx * 4 + 0] * V_port_temp_B[0] + d_G_trans[idx * 4 + 2] * V_port_temp_B[1];
    d_J_trans_B_PVunit_d[idx * 2 + 1] = d_G_trans[idx * 4 + 1] * V_port_temp_B[0] + d_G_trans[idx * 4 + 3] * V_port_temp_B[1];

    //d_J_trans_C_array_d
    double V_port_temp_C[2];
    V_port_temp_C[0] = V_port_C[0] + d_parB_trans[idx * 4 + 0] * I_port_C_PVunit[idx * 2 + 0] + d_parB_trans[idx * 4 + 2] * I_port_C_PVunit[idx * 2 + 1];
    V_port_temp_C[1] = V_port_C[1] + d_parB_trans[idx * 4 + 1] * I_port_C_PVunit[idx * 2 + 0] + d_parB_trans[idx * 4 + 3] * I_port_C_PVunit[idx * 2 + 1];
    d_J_trans_C_PVunit_d[idx * 2 + 0] = d_G_trans[idx * 4 + 0] * V_port_temp_C[0] + d_G_trans[idx * 4 + 2] * V_port_temp_C[1];
    d_J_trans_C_PVunit_d[idx * 2 + 1] = d_G_trans[idx * 4 + 1] * V_port_temp_C[0] + d_G_trans[idx * 4 + 3] * V_port_temp_C[1];
    
	//Capacitor
	for(int i = 0; i < 3; i++)
	{
		C_filt[idx * 3 + i].cap_jh = - C_filt[idx * 3 + i].cap_g * C_filt[idx * 3 + i].cap_v - C_filt[idx * 3 + i].cap_ib;
	}

	//VSC
	double idc_source;
	double jabc[3];

	if(_Dblk1[idx] == 0)
	{
		//DC side
		idc_source = 0;

		//AC side
		jabc[0] = (d_Vin_d[idx][9] - d_Vin_d[idx][10]) * VSC_average[idx].g_vsc;
		jabc[1] = (d_Vin_d[idx][10] - d_Vin_d[idx][11]) * VSC_average[idx].g_vsc;
		jabc[2] = (d_Vin_d[idx][11] - d_Vin_d[idx][9]) * VSC_average[idx].g_vsc;

		//History Current
		for(int i = 0; i < 3; i++)
		{
			L_filt[idx * 3 + i].ind_jh = 0;
		}
	}
	else
	{
		
		idc_source = - SW_VSC[0] * L_filt[idx * 3 + 0].ind_ib - SW_VSC[1] * L_filt[idx * 3 + 1].ind_ib - SW_VSC[2] * L_filt[idx * 3 + 2].ind_ib;

		jabc[0] = (SW_VSC[0] - SW_VSC[1]) * (_VdcBus[idx]) * VSC_average[idx].g_vsc;
		jabc[1] = (SW_VSC[1] - SW_VSC[2]) * (_VdcBus[idx]) * VSC_average[idx].g_vsc;
		jabc[2] = (SW_VSC[2] - SW_VSC[0]) * (_VdcBus[idx]) * VSC_average[idx].g_vsc;

		//Inductor
		for(int i = 0; i < 3; i++)
		{
			L_filt[idx * 3 + i].ind_jh = L_filt[idx * 3 + i].ind_g * L_filt[idx * 3 + i].ind_v + L_filt[idx * 3 + i].ind_ib;
		}
	}

    //OUTPUT
    _idc_source[idx] = idc_source; 

	//Gernate Yg
	Bre_PVunit[idx].g = SW_PVunit * Bre_PVunit[idx].goff + (1 - SW_PVunit) * Bre_PVunit[idx].gon;
	
	d_Yg1_d[idx][18 + 18 * BRANCH_NUM_AC] = Bre_PVunit[idx].g;
	d_Yg1_d[idx][19 + 19 * BRANCH_NUM_AC] = Bre_PVunit[idx].g;
	d_Yg1_d[idx][20 + 20 * BRANCH_NUM_AC] = Bre_PVunit[idx].g;
	
	//Generate J
    for(int i = 0; i < BRANCH_NUM_AC; i++)
    {
        d_J_branch1_d[idx][i] = 0;
    }
	
	d_J_branch1_d[idx][0] = d_J_trans_A_PVunit_d[idx * 2 + 1];
	d_J_branch1_d[idx][5] = d_J_trans_A_PVunit_d[idx * 2 + 0];
	d_J_branch1_d[idx][6] = d_J_trans_B_PVunit_d[idx * 2 + 1];
	d_J_branch1_d[idx][11] = d_J_trans_B_PVunit_d[idx * 2 + 0];
	d_J_branch1_d[idx][12] = d_J_trans_C_PVunit_d[idx * 2 + 1];
	d_J_branch1_d[idx][17] = d_J_trans_C_PVunit_d[idx * 2 + 0];

	d_J_branch1_d[idx][24] = C_filt[idx * 3 + 0].cap_jh;
	d_J_branch1_d[idx][25] = C_filt[idx * 3 + 1].cap_jh;
	d_J_branch1_d[idx][26] = C_filt[idx * 3 + 2].cap_jh;

	d_J_branch1_d[idx][30] = L_filt[idx * 3 + 0].ind_jh;
	d_J_branch1_d[idx][31] = L_filt[idx * 3 + 1].ind_jh;
	d_J_branch1_d[idx][32] = L_filt[idx * 3 + 2].ind_jh;

    d_J_branch1_d[idx][33] = -jabc[0];
	d_J_branch1_d[idx][34] = -jabc[1];
	d_J_branch1_d[idx][35] = -jabc[2];
}

__device__ void PImodel_pre(PImodel* pimodel, double f, double l, double rp, double xp, double xcp, double rz, double xz, double xcz, double delt)
{

    double invdt = 1 / delt;

    double Rp, Xp, Xcp;
    double Rz, Xz, Xcz;

    Rp = rp * l;
    Xp = xp * l;
    Xcp = xcp / l;
    Rz = rz * l;
    Xz = xz * l;
    Xcz = xcz / l;

    //SI units
    double Rs, Rm, Ls, Lm, Cl, Cg, Cm, C;
    Rs = (Rz + 2.0 * Rp) / 3.0;
    Rm = (Rz - Rp) / 3.0;
    Ls = (Xz / (2.0 * PI * f) + 2.0 * Xp / (2.0 * PI * f)) / 3.0;
    Lm = (Xz / (2.0 * PI * f) - Xp / (2.0 * PI * f)) / 3.0;
    Cl = (1.0 / (Xcp * 2.0E6 * PI * f) - 1.0 / (Xcz * 2.0E6 * PI * f)) / 3.0;
    Cg = 1.0 / (Xcz * 2.0E6 * PI * f);

    Cm = Cg / 2.0;
    C = Cl / 2.0;

    double R_series_s = Rs + 2 * Ls * invdt;
    double R_series_m = Rm + 2 * Lm * invdt;
    pimodel->R_cal_s = 2 * Ls * invdt - Rs;
    pimodel->R_cal_m = 2 * Lm * invdt - Rm;

    pimodel->G_series_s = (R_series_s * R_series_s - R_series_m * R_series_m)
        / (R_series_s * R_series_s * R_series_s - 3 * R_series_s * R_series_m * R_series_m + 2 * R_series_m * R_series_m * R_series_m);
    pimodel->G_series_m = R_series_m * (R_series_m - R_series_s)
        / (R_series_s * R_series_s * R_series_s - 3 * R_series_s * R_series_m * R_series_m + 2 * R_series_m * R_series_m * R_series_m);

    //Trapezodial integration of Capacitor
    pimodel->GC = 2.0 * C * invdt;
    pimodel->GCm = 2.0 * Cm * invdt;

    double G1 = pimodel->G_series_s + 2 * pimodel->GC + pimodel->GCm;
    double G2 = pimodel->G_series_m - pimodel->GC;
    double G3 = -pimodel->G_series_s;
    double G4 = -pimodel->G_series_m;

    for(int i = 0; i < 3; i++)
    {
        for(int j = 0; j < 3; j++)
        {
            if(i == j)
            {
                pimodel->G_PImodel[i][j] = G1;
                pimodel->G_PImodel[i + 3][j + 3] = G1;
            }
            else
            {
                pimodel->G_PImodel[i][j] = G2;
                pimodel->G_PImodel[i + 3][j + 3] = G2;
            }
        }
    }

    for(int i = 0; i < 3; i++)
    {
        for(int j = 0; j < 3; j++)
        {
            if(i == j)
            {
                pimodel->G_PImodel[i + 3][j] = G3;
                pimodel->G_PImodel[i][j + 3] = G3;
            }
            else
            {
                pimodel->G_PImodel[i + 3][j] = G4;
                pimodel->G_PImodel[i][j + 3] = G4;
            }
        }
    }
}

__device__ void Trans_pre(double* d_G_trans_sys, double* d_parB_trans_sys, double* Trans_sys_pars, double delt)
{
    double Sb_trans, fb_trans, XL_trans, Z_trans, Im_trans, V1_trans, V2_trans;
    Sb_trans = Trans_sys_pars[0];
	fb_trans = Trans_sys_pars[1];
	XL_trans = Trans_sys_pars[2];
	Z_trans = Trans_sys_pars[3];
	Im_trans = Trans_sys_pars[4]; 
	V1_trans = Trans_sys_pars[5];
	V2_trans = Trans_sys_pars[6];

    double N_trans = V1_trans / V2_trans;
	double Zb1_trans = V1_trans * V1_trans / Sb_trans;
	double Zb2_trans = V2_trans * V2_trans / Sb_trans;

	double L1_trans = (XL_trans / 2) * Zb1_trans / (2 * PI * fb_trans);
	double L2_trans = L1_trans / N_trans / N_trans;
	double Lm_trans = 1.0 / Im_trans * Zb1_trans / (2 * PI * fb_trans); 

	double R1_trans = Z_trans / 2 * Zb1_trans;
	double R2_trans = Z_trans / 2 * Zb2_trans;

    double R_trans[2][2];
    double L_trans[2][2];

    R_trans[0][0] = R1_trans;
    R_trans[0][1] = 0.0;
    R_trans[1][0] = 0.0;
	R_trans[1][1] = R2_trans;

	L_trans[0][0] = L1_trans + Lm_trans;
	L_trans[0][1] = Lm_trans / N_trans;
	L_trans[1][0] = Lm_trans / N_trans;
	L_trans[1][1] = L2_trans + Lm_trans / N_trans / N_trans;

    double parA_trans[2][2];
    for(int i = 0; i < 2; i++)
    {
        for(int j = 0; j < 2; j++)
        {
            parA_trans[i][j] = R_trans[i][j] + 2 / delt * L_trans[i][j];
        }
    }

    double inva = parA_trans[0][0];
    double invb = parA_trans[0][1];
    double invc = parA_trans[1][0];
    double invd = parA_trans[1][1];
    double invdet = inva * invd - invb * invc;
    double inva1 = invd / invdet;
    double invb1 = -invb / invdet;
    double invc1 = -invc / invdet;
    double invd1 = inva / invdet;

    double G_trans[2][2];
    G_trans[0][0] = inva1;
    G_trans[0][1] = invb1;
    G_trans[1][0] = invc1;
    G_trans[1][1] = invd1;

	//ATTENTION 行优先还是列优先
    for (int i = 0; i < 2; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            d_G_trans_sys[i + j * 2] = G_trans[i][j];
        }
    }

    for (int i = 0; i < 2; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            d_parB_trans_sys[i + j * 2] = - R_trans[i][j] + 2 / delt * L_trans[i][j];
        }
    }

}

__global__ void system_pre(double* _DELT, double* _Source_pars, double* _Fault_pars, double* _PImodel_1_pars, double* _Trans_sys_pars, double* _PImodel_2_pars,
    Inductor* L_source, Fault* Fault_1, Fault* Fault_2, PImodel* PImodel_1, PImodel* PImodel_2, double* d_G_trans_sys, double* d_parB_trans_sys)
{
    //Source Inductor
    double L;
    L = _Source_pars[3];

    for (int i = 0; i < 3; i++)
	{
		L_source[i].ind_g = *_DELT / (2.0 * L);
		L_source[i].ind_v = 0.0;
		L_source[i].ind_ib = 0.0;
		L_source[i].ind_jh = 0.0;
	}

    //Fault
    double Ron, Roff;
    Ron = _Fault_pars[0];
    Roff = _Fault_pars[1];
    
    //CHECK!!!
    Fault_1->g_on = 1 / Ron;
    Fault_1->g_off = 1 / Roff;
    Fault_1->g_fault = Fault_1->g_off;

    for(int i = 0; i < 3; i++)
    {
        Fault_2[i].g_on = 1 / Ron;
        Fault_2[i].g_off = 1 / Roff;
        Fault_2[i].g_fault = Fault_2[i].g_off;
    }

    //PImodel
    double f_1, l_1, rp_1, xp_1, xcp_1, rz_1, xz_1, xcz_1;
    double f_2, l_2, rp_2, xp_2, xcp_2, rz_2, xz_2, xcz_2;

    f_1 = _PImodel_1_pars[0];
    l_1 = _PImodel_1_pars[1];
    rp_1 = _PImodel_1_pars[2];
    xp_1 = _PImodel_1_pars[3];
    xcp_1 = _PImodel_1_pars[4];
    rz_1 = _PImodel_1_pars[5];
    xz_1 = _PImodel_1_pars[6];
    xcz_1 = _PImodel_1_pars[7];

    f_2 = _PImodel_2_pars[0];
    l_2 = _PImodel_2_pars[1];
    rp_2 = _PImodel_2_pars[2];
    xp_2 = _PImodel_2_pars[3];
    xcp_2 = _PImodel_2_pars[4];
    rz_2 = _PImodel_2_pars[5];
    xz_2 = _PImodel_2_pars[6];
    xcz_2 = _PImodel_2_pars[7];

    PImodel_pre(PImodel_1, f_1, l_1, rp_1, xp_1, xcp_1, rz_1, xz_1, xcz_1, *_DELT);
    PImodel_pre(PImodel_2, f_2, l_2, rp_2, xp_2, xcp_2, rz_2, xz_2, xcz_2, *_DELT);

    //Transformer
    Trans_pre(d_G_trans_sys, d_parB_trans_sys, _Trans_sys_pars, *_DELT);

}

__global__ void system_emt(double* _Source_pars, double* simulation_time, double* fault1_time, double* fault2_time, Inductor* L_source, Fault* Fault_1, Fault* Fault_2, PImodel* PImodel_1, PImodel* PImodel_2,
    double* d_V_total_d, double* d_Yg_CAL_d, double* d_J_branch_CAL_d, double* G_PVunit, double* INJ_ABC,
    double* d_G_trans_sys, double* d_parB_trans_sys, double* d_J_trans_A_sys_d, double* d_J_trans_B_sys_d, double* d_J_trans_C_sys_d,
    double* d_Ih_mutual1_1, double* d_Ih_mutual2_1, double* d_Ih_self1_1, double* d_Ih_self2_1, double* d_Ih_ABC_1,
    double* d_Ih_mutual1_2, double* d_Ih_mutual2_2, double* d_Ih_self1_2, double* d_Ih_self2_2, double* d_Ih_ABC_2,
    double* d_Pstation, double* output_his_LPF_Pstation, double* input_his_LPF_Pstation, double* d_Qstation, double* output_his_LPF_Qstation, double* input_his_LPF_Qstation, double* d_Vstation, double* d_Istation)
{
    //Source
    double Vm, freq, Rin;
    Vm = _Source_pars[0];
    freq = _Source_pars[1];
    Rin = _Source_pars[2];

    double G_source = 1 / Rin;

    double Mag = Vm * sqrt_2 / sqrt_3;
    double w = 2 * PI * freq;

    double phaseShiftInRadians[3];
    phaseShiftInRadians[0] = 0 / 180.0 * PI;
    phaseShiftInRadians[1] = -120.0 / 180.0 * PI;
    phaseShiftInRadians[2] = 120.0 / 180.0 * PI;

    double J_source[3];
    for(int i = 0; i < 3; i++)
    {
        J_source[i] = Mag * sin(w * (*simulation_time) + phaseShiftInRadians[i]) * G_source;
    }

    //Inductor     
	L_source[0].ind_v = d_V_total_d[0] - d_V_total_d[3];
	L_source[1].ind_v = d_V_total_d[1] - d_V_total_d[4];
	L_source[2].ind_v = d_V_total_d[2] - d_V_total_d[5];

    for(int i = 0; i < 3; i++)
	{
		L_source[i].ind_ib = L_source[i].ind_g * L_source[i].ind_v + L_source[i].ind_jh;
	}

    for(int i = 0; i < 3; i++)
	{
		L_source[i].ind_jh = L_source[i].ind_g * L_source[i].ind_v + L_source[i].ind_ib;
	}

    //Fault1 & Fault2
    if(*simulation_time >= *fault1_time && *simulation_time <= *fault1_time + 0.10346)
    {
        Fault_1->g_fault = Fault_1->g_on;
    }
    else
    {
        Fault_1->g_fault = Fault_1->g_off;
    }

    if(*simulation_time >= *fault2_time && *simulation_time <= *fault2_time + 0.1)
    {
        Fault_2[0].g_fault = Fault_2[0].g_on;
        Fault_2[1].g_fault = Fault_2[1].g_on;
        Fault_2[2].g_fault = Fault_2[2].g_on;
    }
    else
    {
        Fault_2[0].g_fault = Fault_2[0].g_off;
        Fault_2[1].g_fault = Fault_2[1].g_off;
        Fault_2[2].g_fault = Fault_2[2].g_off;
    }

    //PImodel_1
    double Vabc1_1[3], Vabc2_1[3], V_ABC_1[3];
        
    Vabc1_1[0] = d_V_total_d[3] - d_V_total_d[4];
    Vabc1_1[1] = d_V_total_d[4] - d_V_total_d[5];
    Vabc1_1[2] = d_V_total_d[5] - d_V_total_d[3];
        
    Vabc2_1[0] = d_V_total_d[6] - d_V_total_d[7];
    Vabc2_1[1] = d_V_total_d[7] - d_V_total_d[8];
    Vabc2_1[2] = d_V_total_d[8] - d_V_total_d[6];

    for(int i = 0; i < 3; i++)
    {
        V_ABC_1[i] = d_V_total_d[i + 3] - d_V_total_d[i + 6];
    }

    double Imutual1_1[3], Imutual2_1[3], Iself1_1[3], Iself2_1[3], I_ABC_1[3];
    for(int i = 0; i < 3; i ++)
    {
        Imutual1_1[i] = Vabc1_1[i] * PImodel_1->GC + d_Ih_mutual1_1[i];
        Imutual2_1[i] = Vabc2_1[i] * PImodel_1->GC + d_Ih_mutual2_1[i];

        Iself1_1[i] = d_V_total_d[i + 3] * PImodel_1->GCm + d_Ih_self1_1[i];
        Iself2_1[i] = d_V_total_d[i + 6] * PImodel_1->GCm + d_Ih_self2_1[i];
    }

    I_ABC_1[0] = PImodel_1->G_series_s * V_ABC_1[0] + PImodel_1->G_series_m * V_ABC_1[1] + PImodel_1->G_series_m * V_ABC_1[2] + d_Ih_ABC_1[0];
    I_ABC_1[1] = PImodel_1->G_series_m * V_ABC_1[0] + PImodel_1->G_series_s * V_ABC_1[1] + PImodel_1->G_series_m * V_ABC_1[2] + d_Ih_ABC_1[1];
    I_ABC_1[2] = PImodel_1->G_series_m * V_ABC_1[0] + PImodel_1->G_series_m * V_ABC_1[1] + PImodel_1->G_series_s * V_ABC_1[2] + d_Ih_ABC_1[2];

    for(int i = 0; i < 3; i ++)
    {
        d_Ih_mutual1_1[i] = - Vabc1_1[i] * PImodel_1->GC - Imutual1_1[i];
        d_Ih_mutual2_1[i] = - Vabc2_1[i] * PImodel_1->GC - Imutual2_1[i];

        d_Ih_self1_1[i] = - d_V_total_d[i + 3] * PImodel_1->GCm - Iself1_1[i];
        d_Ih_self2_1[i] = - d_V_total_d[i + 6] * PImodel_1->GCm - Iself2_1[i];
    }

    double UI_a_1, UI_b_1, UI_c_1;
    UI_a_1 = V_ABC_1[0] + PImodel_1->R_cal_s * I_ABC_1[0] + PImodel_1->R_cal_m * I_ABC_1[1] + PImodel_1->R_cal_m * I_ABC_1[2];
    UI_b_1 = V_ABC_1[1] + PImodel_1->R_cal_m * I_ABC_1[0] + PImodel_1->R_cal_s * I_ABC_1[1] + PImodel_1->R_cal_m * I_ABC_1[2];
    UI_c_1 = V_ABC_1[2] + PImodel_1->R_cal_m * I_ABC_1[0] + PImodel_1->R_cal_m * I_ABC_1[1] + PImodel_1->R_cal_s * I_ABC_1[2];

    d_Ih_ABC_1[0] = PImodel_1->G_series_s * UI_a_1 + PImodel_1->G_series_m * UI_b_1 + PImodel_1->G_series_m * UI_c_1;
    d_Ih_ABC_1[1] = PImodel_1->G_series_m * UI_a_1 + PImodel_1->G_series_s * UI_b_1 + PImodel_1->G_series_m * UI_c_1;
    d_Ih_ABC_1[2] = PImodel_1->G_series_m * UI_a_1 + PImodel_1->G_series_m * UI_b_1 + PImodel_1->G_series_s * UI_c_1;

    //Transformer
	double V_port_A[2], V_port_B[2], V_port_C[2];
    double I_port_A[2], I_port_B[2], I_port_C[2];

	V_port_A[0] = d_V_total_d[6] - 0;
	V_port_B[0] = d_V_total_d[7] - 0;
	V_port_C[0] = d_V_total_d[8] - 0;

	V_port_A[1] = d_V_total_d[9] - d_V_total_d[10];
	V_port_B[1] = d_V_total_d[10] - d_V_total_d[11];
	V_port_C[1] = d_V_total_d[11] - d_V_total_d[9];

    I_port_A[0] = d_G_trans_sys[0] * V_port_A[0] + d_G_trans_sys[2] * V_port_A[1] + d_J_trans_A_sys_d[0];
    I_port_A[1] = d_G_trans_sys[1] * V_port_A[0] + d_G_trans_sys[3] * V_port_A[1] + d_J_trans_A_sys_d[1];
	I_port_B[0] = d_G_trans_sys[0] * V_port_B[0] + d_G_trans_sys[2] * V_port_B[1] + d_J_trans_B_sys_d[0];
    I_port_B[1] = d_G_trans_sys[1] * V_port_B[0] + d_G_trans_sys[3] * V_port_B[1] + d_J_trans_B_sys_d[1];
    I_port_C[0] = d_G_trans_sys[0] * V_port_C[0] + d_G_trans_sys[2] * V_port_C[1] + d_J_trans_C_sys_d[0];
    I_port_C[1] = d_G_trans_sys[1] * V_port_C[0] + d_G_trans_sys[3] * V_port_C[1] + d_J_trans_C_sys_d[1];

    //P and Q
    double Kf = 2 * 0.01 / SIMULATION_DELT;
    double Pstation, Qstation;

    d_Vstation[0] = d_V_total_d[9];
    d_Vstation[1] = d_V_total_d[10];
    d_Vstation[2] = d_V_total_d[11];

    d_Istation[0] = I_port_A[1] - I_port_C[1];
    d_Istation[1] = I_port_B[1] - I_port_A[1];
    d_Istation[2] = I_port_C[1] - I_port_B[1];

    Pstation = d_Vstation[0] * d_Istation[0] + d_Vstation[1] * d_Istation[1] + d_Vstation[2] * d_Istation[2];
    LPF_Execute(Kf, 1.0, Pstation, output_his_LPF_Pstation, input_his_LPF_Pstation, d_Pstation);

    Qstation = ((d_Vstation[1] - d_Vstation[2]) * d_Istation[0] + (d_Vstation[2] - d_Vstation[0]) * d_Istation[1] + (d_Vstation[0] - d_Vstation[1]) * d_Istation[2]) / 1.7320508;
    LPF_Execute(Kf, 1.0, Qstation, output_his_LPF_Qstation, input_his_LPF_Qstation, d_Qstation);

    //d_J_trans_A_sys_d
    double V_port_temp_A[2];
    V_port_temp_A[0] = V_port_A[0] + d_parB_trans_sys[0] * I_port_A[0] + d_parB_trans_sys[2] * I_port_A[1];
    V_port_temp_A[1] = V_port_A[1] + d_parB_trans_sys[1] * I_port_A[0] + d_parB_trans_sys[3] * I_port_A[1];
    d_J_trans_A_sys_d[0] = d_G_trans_sys[0] * V_port_temp_A[0] + d_G_trans_sys[2] * V_port_temp_A[1];
    d_J_trans_A_sys_d[1] = d_G_trans_sys[1] * V_port_temp_A[0] + d_G_trans_sys[3] * V_port_temp_A[1];

    //d_J_trans_B_sys_d
    double V_port_temp_B[2];
    V_port_temp_B[0] = V_port_B[0] + d_parB_trans_sys[0] * I_port_B[0] + d_parB_trans_sys[2] * I_port_B[1];
    V_port_temp_B[1] = V_port_B[1] + d_parB_trans_sys[1] * I_port_B[0] + d_parB_trans_sys[3] * I_port_B[1];
    d_J_trans_B_sys_d[0] = d_G_trans_sys[0] * V_port_temp_B[0] + d_G_trans_sys[2] * V_port_temp_B[1];
    d_J_trans_B_sys_d[1] = d_G_trans_sys[1] * V_port_temp_B[0] + d_G_trans_sys[3] * V_port_temp_B[1];

    //d_J_trans_C_sys_d
    double V_port_temp_C[2];
    V_port_temp_C[0] = V_port_C[0] + d_parB_trans_sys[0] * I_port_C[0] + d_parB_trans_sys[2] * I_port_C[1];
    V_port_temp_C[1] = V_port_C[1] + d_parB_trans_sys[1] * I_port_C[0] + d_parB_trans_sys[3] * I_port_C[1];
    d_J_trans_C_sys_d[0] = d_G_trans_sys[0] * V_port_temp_C[0] + d_G_trans_sys[2] * V_port_temp_C[1];
    d_J_trans_C_sys_d[1] = d_G_trans_sys[1] * V_port_temp_C[0] + d_G_trans_sys[3] * V_port_temp_C[1];

    //PImodel_2
    double Vabc1_2[4][3], Vabc2_2[4][3], V_ABC_2[4][3];
        
    //i represents the number of farm
    for(int i = 0; i < 4; i++)
    {
        Vabc1_2[i][0] = d_V_total_d[9] - d_V_total_d[10];
        Vabc1_2[i][1] = d_V_total_d[10] - d_V_total_d[11];
        Vabc1_2[i][2] = d_V_total_d[11] - d_V_total_d[9];
        
        Vabc2_2[i][0] = d_V_total_d[i * 3 + 12] - d_V_total_d[i * 3 + 13];
        Vabc2_2[i][1] = d_V_total_d[i * 3 + 13] - d_V_total_d[i * 3 + 14];
        Vabc2_2[i][2] = d_V_total_d[i * 3 + 14] - d_V_total_d[i * 3 + 12];
    }
        
    for(int i = 0; i < 4; i++)
    {
        for(int j = 0; j < 3; j++)
        {
            V_ABC_2[i][j] = d_V_total_d[j + 9]- d_V_total_d[i * 3 + j + 12];
        }    
    }

    double Imutual1_2[4][3], Imutual2_2[4][3], Iself1_2[4][3], Iself2_2[4][3], I_ABC_2[4][3];
        
    for(int i = 0; i < 4; i ++)
    {
        for(int j = 0; j < 3; j++)
        {
            Imutual1_2[i][j] = Vabc1_2[i][j] * PImodel_2->GC + d_Ih_mutual1_2[i * 4 + j];
            Imutual2_2[i][j] = Vabc2_2[i][j] * PImodel_2->GC + d_Ih_mutual2_2[i * 4 + j];

            Iself1_2[i][j] = d_V_total_d[j + 9] * PImodel_2->GCm + d_Ih_self1_2[i * 4 + j];
            Iself2_2[i][j] = d_V_total_d[i * 3 + j + 12] * PImodel_2->GCm + d_Ih_self2_2[i * 4 + j];
        }  
    }

    for(int i = 0; i < 4; i++)
    {
        I_ABC_2[i][0] = PImodel_2->G_series_s * V_ABC_2[i][0] + PImodel_2->G_series_m * V_ABC_2[i][1] + PImodel_2->G_series_m * V_ABC_2[i][2] + d_Ih_ABC_2[i * 4 + 0];
        I_ABC_2[i][1] = PImodel_2->G_series_m * V_ABC_2[i][0] + PImodel_2->G_series_s * V_ABC_2[i][1] + PImodel_2->G_series_m * V_ABC_2[i][2] + d_Ih_ABC_2[i * 4 + 1];
        I_ABC_2[i][2] = PImodel_2->G_series_m * V_ABC_2[i][0] + PImodel_2->G_series_m * V_ABC_2[i][1] + PImodel_2->G_series_s * V_ABC_2[i][2] + d_Ih_ABC_2[i * 4 + 2];
    }
        
    for(int i = 0; i < 4; i++)
    {
        for(int j = 0; j < 3; j ++)
        {
            d_Ih_mutual1_2[i * 4 + j] = - Vabc1_2[i][j] * PImodel_2->GC - Imutual1_2[i][j];
            d_Ih_mutual2_2[i * 4 + j] = - Vabc2_2[i][j] * PImodel_2->GC - Imutual2_2[i][j];

            d_Ih_self1_2[i * 4 + j] = - d_V_total_d[j + 9] * PImodel_2->GCm - Iself1_2[i][j];
            d_Ih_self2_2[i * 4 + j] = - d_V_total_d[i * 3 + j + 12] * PImodel_2->GCm - Iself2_2[i][j];
        }
    }

    double UI_a_2[4], UI_b_2[4], UI_c_2[4];
        
    for(int i = 0; i < 4; i++)
    {
        UI_a_2[i] = V_ABC_2[i][0] + PImodel_2->R_cal_s * I_ABC_2[i][0] + PImodel_2->R_cal_m * I_ABC_2[i][1] + PImodel_2->R_cal_m * I_ABC_2[i][2];
        UI_b_2[i] = V_ABC_2[i][1] + PImodel_2->R_cal_m * I_ABC_2[i][0] + PImodel_2->R_cal_s * I_ABC_2[i][1] + PImodel_2->R_cal_m * I_ABC_2[i][2];
        UI_c_2[i] = V_ABC_2[i][2] + PImodel_2->R_cal_m * I_ABC_2[i][0] + PImodel_2->R_cal_m * I_ABC_2[i][1] + PImodel_2->R_cal_s * I_ABC_2[i][2];

        d_Ih_ABC_2[i * 4 + 0] = PImodel_2->G_series_s * UI_a_2[i] + PImodel_2->G_series_m * UI_b_2[i] + PImodel_2->G_series_m * UI_c_2[i];
        d_Ih_ABC_2[i * 4 + 1] = PImodel_2->G_series_m * UI_a_2[i] + PImodel_2->G_series_s * UI_b_2[i] + PImodel_2->G_series_m * UI_c_2[i];
        d_Ih_ABC_2[i * 4 + 2] = PImodel_2->G_series_m * UI_a_2[i] + PImodel_2->G_series_m * UI_b_2[i] + PImodel_2->G_series_s * UI_c_2[i];
    }

    /////////////////////////Generate Yg/////////////////////////////////////
	d_Yg_CAL_d[0 + 0 * BRANCH_NUM] = G_source;
	d_Yg_CAL_d[1 + 1 * BRANCH_NUM] = G_source;
	d_Yg_CAL_d[2 + 2 * BRANCH_NUM] = G_source;

	d_Yg_CAL_d[3 + 3 * BRANCH_NUM] = L_source[0].ind_g;
	d_Yg_CAL_d[4 + 4 * BRANCH_NUM] = L_source[1].ind_g;
	d_Yg_CAL_d[5 + 5 * BRANCH_NUM] = L_source[2].ind_g;

	//three-phase fault
    d_Yg_CAL_d[6 + 6 * BRANCH_NUM] = Fault_1->g_fault;
	d_Yg_CAL_d[7 + 7 * BRANCH_NUM] = Fault_1->g_off;
	d_Yg_CAL_d[8 + 8 * BRANCH_NUM] = Fault_1->g_off;

    //Transformer
	for(int i = 0; i < 3; i++)
	{
		d_Yg_CAL_d[i * 6 + 9 + (i * 6 + 9) * BRANCH_NUM] = d_G_trans_sys[0 + 0 * 2];
		d_Yg_CAL_d[i * 6 + 10 + (i * 6 + 10) * BRANCH_NUM] = -d_G_trans_sys[0 + 1 * 2];
		d_Yg_CAL_d[i * 6 + 11 + (i * 6 + 11) * BRANCH_NUM] = d_G_trans_sys[0 + 1 * 2];
		d_Yg_CAL_d[i * 6 + 12 + (i * 6 + 12) * BRANCH_NUM] = -d_G_trans_sys[0 + 1 * 2];
		d_Yg_CAL_d[i * 6 + 13 + (i * 6 + 13) * BRANCH_NUM] = d_G_trans_sys[0 + 1 * 2];
		d_Yg_CAL_d[i * 6 + 14 + (i * 6 + 14) * BRANCH_NUM] = d_G_trans_sys[1 + 1 * 2];
	}

    //single-phase fault
    d_Yg_CAL_d[27 + 27 * BRANCH_NUM] = Fault_2[0].g_fault;
    d_Yg_CAL_d[28 + 28 * BRANCH_NUM] = Fault_2[1].g_fault;
    d_Yg_CAL_d[29 + 29 * BRANCH_NUM] = Fault_2[2].g_fault;

    //Total G of PVFarm
    double G_PVFarm[4][6] = {0};

    for(int i = 0; i < 6; i++)
    {
        for(int j = 0; j < DIM; j++)
        {
            G_PVFarm[0][i] += G_PVunit[j * 6 + i];
            G_PVFarm[1][i] += G_PVunit[(j + DIM) * 6 + i];
            G_PVFarm[2][i] += G_PVunit[(j + DIM * 2) * 6 + i];
            G_PVFarm[3][i] += G_PVunit[(j + DIM * 3) * 6 + i];
        }
    }
        
    for(int i = 0; i < 4; i++)
    {
        d_Yg_CAL_d[i * 6 + 30 + (i * 6 + 30) * BRANCH_NUM] = G_PVFarm[i][3];
        d_Yg_CAL_d[i * 6 + 31 + (i * 6 + 31) * BRANCH_NUM] = G_PVFarm[i][4];
        d_Yg_CAL_d[i * 6 + 32 + (i * 6 + 32) * BRANCH_NUM] = G_PVFarm[i][5];

        d_Yg_CAL_d[i * 6 + 33 + (i * 6 + 33) * BRANCH_NUM] = G_PVFarm[i][0];
        d_Yg_CAL_d[i * 6 + 34 + (i * 6 + 34) * BRANCH_NUM] = G_PVFarm[i][1];
        d_Yg_CAL_d[i * 6 + 35 + (i * 6 + 35) * BRANCH_NUM] = G_PVFarm[i][2];
    }

    /////////////////////////Generate J////////////////////////////
    d_J_branch_CAL_d[0] = J_source[0];
    d_J_branch_CAL_d[1] = J_source[1];
    d_J_branch_CAL_d[2] = J_source[2];
        
    d_J_branch_CAL_d[3]= L_source[0].ind_jh;
    d_J_branch_CAL_d[4] = L_source[1].ind_jh;
    d_J_branch_CAL_d[5] = L_source[2].ind_jh;
        
    d_J_branch_CAL_d[9] = d_J_trans_A_sys_d[0];
	d_J_branch_CAL_d[14] = d_J_trans_A_sys_d[1];
	d_J_branch_CAL_d[15] = d_J_trans_B_sys_d[0];
	d_J_branch_CAL_d[20] = d_J_trans_B_sys_d[1];
    d_J_branch_CAL_d[21] = d_J_trans_C_sys_d[0];
    d_J_branch_CAL_d[26] = d_J_trans_C_sys_d[1];

    //Total INJ of PVFarm
    double INJ_PVFarm[4][3] = {0};

    for(int i = 0; i < 3; i++)
    {
        for(int j = 0; j < DIM; j++)
        {
            INJ_PVFarm[0][i] += INJ_ABC[j * 3 + i];
            INJ_PVFarm[1][i] += INJ_ABC[(j + DIM) * 3 + i];
            INJ_PVFarm[2][i] += INJ_ABC[(j + DIM * 2) * 3 + i];
            INJ_PVFarm[3][i] += INJ_ABC[(j + DIM * 3) * 3 + i];
        }
    }
        
    for(int i = 0; i < 4; i++)
    {
        d_J_branch_CAL_d[i * 6 + 33] = INJ_PVFarm[i][0];
        d_J_branch_CAL_d[i * 6 + 34] = INJ_PVFarm[i][1];
        d_J_branch_CAL_d[i * 6 + 35] = INJ_PVFarm[i][2];
    }
            
}

__global__ void Y_total_Calculation(double* d_Y_total_d, double* d_Y_total_inv_d, double* d_Y_CAL_d, PImodel* PImodel_1, PImodel* PImodel_2)
{
    for(int i = 0; i < NODE_NUM; i++)
    {
        for(int j = 0; j < NODE_NUM; j++)
        {
            d_Y_total_d[i + j * NODE_NUM] = d_Y_CAL_d[i + j * NODE_NUM]; 
        }
    }   

    //PImodel_1
    for(int i = 0; i < 6; i++)
    {
        for(int j = 0; j < 6; j++)
        {
            d_Y_total_d[i + 3 + (j + 3) * NODE_NUM] += PImodel_1->G_PImodel[i][j];
        }
    }

    //PImodel_2
    for(int i = 0; i < 4; i++)
    {
        for(int j = 0; j < 3; j++)
        {
            for(int k = 0; k < 3; k++)
            {
                d_Y_total_d[j + 9 + (k + 9) * NODE_NUM] += PImodel_2->G_PImodel[j][k];
                d_Y_total_d[j + 9 + (i * 3 + k + 12) * NODE_NUM] += PImodel_2->G_PImodel[j][k + 3];
                d_Y_total_d[i * 3 + j + 12 + (k + 9) * NODE_NUM] += PImodel_2->G_PImodel[j + 3][k];
                d_Y_total_d[i * 3 + j + 12 + (i * 3 + k + 12) * NODE_NUM] += PImodel_2->G_PImodel[j + 3][k + 3];
            }
        }
    }

    //Y_total_inv 先单位化
    for(int i = 0; i < NODE_NUM; i++)
    {
        for(int j = 0; j < NODE_NUM; j++)
        {
            d_Y_total_inv_d[i + j * NODE_NUM] = 0.0; 
        }
        
    }  

    for(int i = 0; i < NODE_NUM; i++)
    {
        d_Y_total_inv_d[i + i * NODE_NUM] = 1.0; 
    }   

}

__global__ void J_total_Calculation(double* d_J_total_d, double* d_J_CAL_d,
    double* d_Ih_mutual1_1, double* d_Ih_mutual2_1, double* d_Ih_self1_1, double* d_Ih_self2_1, double* d_Ih_ABC_1,
    double* d_Ih_mutual1_2, double* d_Ih_mutual2_2, double* d_Ih_self1_2, double* d_Ih_self2_2, double* d_Ih_ABC_2)
{
    
    for(int i = 0; i < NODE_NUM; i++)
	{
		d_J_total_d[i] = d_J_CAL_d[i];
	}
    
    //PImodel_1
    d_J_total_d[3] += -d_Ih_mutual1_1[0] + d_Ih_mutual1_1[2] - d_Ih_self1_1[0] - d_Ih_ABC_1[0];
    d_J_total_d[4] += -d_Ih_mutual1_1[1] + d_Ih_mutual1_1[0] - d_Ih_self1_1[1] - d_Ih_ABC_1[1];
    d_J_total_d[5] += -d_Ih_mutual1_1[2] + d_Ih_mutual1_1[1] - d_Ih_self1_1[2] - d_Ih_ABC_1[2];

    d_J_total_d[6] += -d_Ih_mutual2_1[0] + d_Ih_mutual2_1[2] - d_Ih_self2_1[0] + d_Ih_ABC_1[0];
    d_J_total_d[7] += -d_Ih_mutual2_1[1] + d_Ih_mutual2_1[0] - d_Ih_self2_1[1] + d_Ih_ABC_1[1];
    d_J_total_d[8] += -d_Ih_mutual2_1[2] + d_Ih_mutual2_1[1] - d_Ih_self2_1[2] + d_Ih_ABC_1[2];

    //PImodel_2
    for(int i = 0; i < 4; i++)
    {
        d_J_total_d[9]  += -d_Ih_mutual1_2[i * 4 + 0] + d_Ih_mutual1_2[i * 4 + 2] - d_Ih_self1_2[i * 4 + 0] - d_Ih_ABC_2[i * 4 + 0];
        d_J_total_d[10] += -d_Ih_mutual1_2[i * 4 + 1] + d_Ih_mutual1_2[i * 4 + 0] - d_Ih_self1_2[i * 4 + 1] - d_Ih_ABC_2[i * 4 + 1];
        d_J_total_d[11] += -d_Ih_mutual1_2[i * 4 + 2] + d_Ih_mutual1_2[i * 4 + 1] - d_Ih_self1_2[i * 4 + 2] - d_Ih_ABC_2[i * 4 + 2];

        d_J_total_d[i * 3 + 12] += -d_Ih_mutual2_2[i * 4 + 0] + d_Ih_mutual2_2[i * 4 + 2] -d_Ih_self2_2[i * 4 + 0] + d_Ih_ABC_2[i * 4 + 0];
        d_J_total_d[i * 3 + 13] += -d_Ih_mutual2_2[i * 4 + 1] + d_Ih_mutual2_2[i * 4 + 0] -d_Ih_self2_2[i * 4 + 1] + d_Ih_ABC_2[i * 4 + 1];
        d_J_total_d[i * 3 + 14] += -d_Ih_mutual2_2[i * 4 + 2] + d_Ih_mutual2_2[i * 4 + 1] -d_Ih_self2_2[i * 4 + 2] + d_Ih_ABC_2[i * 4 + 2];
    }
}

__global__ void Output_System(double* V_total, double* V_ABC, double** Vex)
{
    int idx = threadIdx.x;

    for(int j = 0; j < 3; j++)
    {
        V_ABC[idx * 3 + j] = V_total[12 + j];
        V_ABC[(idx + DIM) * 3 + j] = V_total[15 + j];
        V_ABC[(idx + DIM * 2) * 3 + j] = V_total[18 + j];
        V_ABC[(idx + DIM * 3) * 3 + j] = V_total[21 + j];
    }

    
    for(int i = 0; i < 4; i++)
    {
        Vex[idx + DIM * i][0] = V_ABC[(idx + DIM * i) * 3 + 0];
        Vex[idx + DIM * i][1] = V_ABC[(idx + DIM * i) * 3 + 1];
        Vex[idx + DIM * i][2] = V_ABC[(idx + DIM * i) * 3 + 2];
    }
    
}

__global__ void Output_DCside(double* V_PV_kV, double* I_PV_kA, double** V_DC, double* VpvBus, double* Ipv, double* VdcBus)
{
    int idx = threadIdx.x;
    
    VpvBus[idx] = V_PV_kV[idx];
    Ipv[idx] = I_PV_kA[idx];
    VdcBus[idx] = V_DC[idx][0] - 0;
    
}

__global__ void Output_ACside(double** Vin, double* I_port_A, double* I_port_B, double* I_port_C,
    double* Vcon1A, double* Vcon1B, double* Vcon1C, double* Icon1A, double* Icon1B, double* Icon1C, double* Icon1A1, double* Icon1B1, double* Icon1C1, 
    double* output_his_LPF_Pcon1, double* input_his_LPF_Pcon1, double* output_his_LPF_Qcon1, double* input_his_LPF_Qcon1, double* output_his_LPF_VacRMS, double* input_his_LPF_VacRMS,
    double* Pcon1filt, double* Qcon1filt, double* VacRMSfilt)
{
    int idx = threadIdx.x;
       
    Vcon1A[idx] = Vin[idx][0];
    Vcon1B[idx] = Vin[idx][1];
    Vcon1C[idx] = Vin[idx][2];
    Icon1A[idx] = I_port_A[idx * 2 + 0];
    Icon1B[idx] = I_port_B[idx * 2 + 0];
    Icon1C[idx] = I_port_C[idx * 2 + 0];
    Icon1A1[idx] = - Icon1A[idx];
    Icon1B1[idx] = - Icon1B[idx];
    Icon1C1[idx] = - Icon1C[idx];

    //LPF
    double Kf = 2 * 0.01 / SIMULATION_DELT;
    double Pcon1, Qcon1, VacRMS;

    Pcon1 = Vcon1A[idx] * Icon1A[idx] + Vcon1B[idx] * Icon1B[idx] + Vcon1C[idx] * Icon1C[idx];
    LPF_Execute(Kf, 1.0, Pcon1, &output_his_LPF_Pcon1[idx], &input_his_LPF_Pcon1[idx], &Pcon1filt[idx]);

    Qcon1 = ((Vcon1B[idx] - Vcon1C[idx]) * Icon1A[idx] + (Vcon1C[idx] - Vcon1A[idx]) * Icon1B[idx] + (Vcon1A[idx] - Vcon1B[idx]) * Icon1C[idx]) / 1.7320508;
    LPF_Execute(Kf, 1.0, Qcon1, &output_his_LPF_Qcon1[idx], &input_his_LPF_Qcon1[idx], &Qcon1filt[idx]);

    VacRMS = sqrt(( pow((Vcon1A[idx] - Vcon1B[idx]) , 2) + pow((Vcon1B[idx] - Vcon1C[idx]) , 2) + pow((Vcon1C[idx] - Vcon1A[idx]) , 2)) / 3);
    LPF_Execute(Kf, 1.0, VacRMS, &output_his_LPF_VacRMS[idx], &input_his_LPF_VacRMS[idx], &VacRMSfilt[idx]);
    
}

__global__ void Output_print(double* simulation_time, int* simulation_num, double* VpvBus, double* Ipv, double* VdcBus, double* Pcon1filt, double* Qcon1filt, double* VacRMSfilt, int *SW_AC, int *SW_DC, double *Pstation, double *Qstation, double *Vstation, double *Istation,
    double* out1, double* out2, double* out3, double* out4, double* out5, double* out6, double* out7,
    double* out_1, double* out_2, double* out_3, double* out_4, double* out_5, double* out_6, double* out_7,
    double* out_control_1, int* out_control_2, int* out_control_3, int* out_control_4, int* out_control_5, int* out_control_6, int* out_control_7, int* out_control_8,
    double* out_system_1, double* out_system_2, double* out_system_3, double* out_system_4, double* out_system_5, double* out_system_6, double* out_system_7, double* out_system_8, double* out_system_9)
{
    
    out1[*simulation_num] = *simulation_time;
    out2[*simulation_num] = VpvBus[0];
    out3[*simulation_num] = Ipv[0];
    out4[*simulation_num] = VdcBus[0];
    out5[*simulation_num] = Pcon1filt[0];
    out6[*simulation_num] = Qcon1filt[0];
    out7[*simulation_num] = VacRMSfilt[0];

    out_1[*simulation_num] = *simulation_time;
    out_2[*simulation_num] = VpvBus[31];
    out_3[*simulation_num] = Ipv[31];
    out_4[*simulation_num] = VdcBus[31];
    out_5[*simulation_num] = Pcon1filt[31];
    out_6[*simulation_num] = Qcon1filt[31];
    out_7[*simulation_num] = VacRMSfilt[31];

    out_control_1[*simulation_num] = *simulation_time;
    out_control_2[*simulation_num] = SW_DC[0];
    out_control_3[*simulation_num] = SW_DC[1];
    out_control_4[*simulation_num] = SW_DC[2];
    out_control_5[*simulation_num] = SW_DC[3];
    out_control_6[*simulation_num] = SW_AC[0];
    out_control_7[*simulation_num] = SW_AC[1];
    out_control_8[*simulation_num] = SW_AC[2];

    out_system_1[*simulation_num] = *simulation_time;
    out_system_2[*simulation_num] = *Pstation;
    out_system_3[*simulation_num] = *Qstation;
    out_system_4[*simulation_num] = Vstation[0];
    out_system_5[*simulation_num] = Vstation[1];
    out_system_6[*simulation_num] = Vstation[2];
    out_system_7[*simulation_num] = Istation[0];
    out_system_8[*simulation_num] = Istation[1];
    out_system_9[*simulation_num] = Istation[2];

    *simulation_time += SIMULATION_DELT;
    *simulation_num += 1;

}

__device__ void exchangeRows(double* mat, int row1, int row2, int n) {
    int tx = threadIdx.x;
    for (int col = tx; col < n; col += blockDim.x) {
        double temp = mat[row1 * n + col];
        mat[row1 * n + col] = mat[row2 * n + col];
        mat[row2 * n + col] = temp;
    }
}

__global__ void invertMatrixKernel(double* d_A, double* d_Ainv, int n) {
    extern __shared__ double shared_data[];

    double* mat = shared_data;       
    double* inv = &shared_data[n * n]; 

    int tx = threadIdx.x;

    for (int i = tx; i < n * n; i += blockDim.x) {
        mat[i] = d_A[i];
        inv[i] = (i / n == i % n) ? 1.0 : 0.0; 
    }
    __syncthreads();

    for (int i = 0; i < n; ++i) {
        if (tx == 0) {
            double max_val = fabs(mat[i * n + i]);
            int max_row = i;
            for (int k = i + 1; k < n; ++k) {
                if (fabs(mat[k * n + i]) > max_val) 
                {
                    max_val = fabs(mat[k * n + i]);
                    max_row = k;
                }
            }

            if (max_row != i) {
                exchangeRows(mat, i, max_row, n);
                exchangeRows(inv, i, max_row, n);
            }
        }
        __syncthreads();


        double diag_val = mat[i * n + i];
        if (diag_val == 0) return; 
        for (int j = tx; j < n; j += blockDim.x) {
            mat[i * n + j] /= diag_val;
            inv[i * n + j] /= diag_val;
        }
        __syncthreads();

        for (int k = 0; k < n; ++k) {
            if (k == i) continue;
            double factor = mat[k * n + i];
            for (int j = tx; j < n; j += blockDim.x) {
                mat[k * n + j] -= factor * mat[i * n + j];
                inv[k * n + j] -= factor * inv[i * n + j];
            }
        }
        __syncthreads();
    }

    for (int i = tx; i < n * n; i += blockDim.x) {
        d_Ainv[i] = inv[i];
    }
}


int main()
{
    ///////////////////////System Parameter//////////////////////////////
    //Source
    double Vm, freq, Rin, L;
    Vm = 220.0; //[kV]
    freq = 50.0; //[Hz]
    Rin = 0.65; //[ohm]
    L = 0.07; //[H]

    //Fault [ABC->G]
    double Ron, Roff;
    Ron = 0.01; //[ohm]
    Roff = 1.0E6; //[ohm]

    //PImodel_1
    double f_1, l_1, rp_1, xp_1, xcp_1, rz_1, xz_1, xcz_1;
    f_1 = 50.0; //[Hz]
    l_1 = 1.0E5; //[m]
    rp_1 = 0.1781598E-4; //[ohm/m]
    xp_1 = 0.31388E-3; //[ohn/m]
    xcp_1 = 273.5448; //[Mohm*m]
    rz_1 = 0.2952201E-3; //[ohn/m]
    xz_1 = 0.1039898E-2; //[ohm/m]
    xcz_1 = 414.1642; //[Mohm*m]

    //Transformer
    double Sb_trans, fb_trans, XL_trans, Z_trans, Im_trans, V1_trans, V2_trans;
    Sb_trans = 33.3; //[MVA]
	fb_trans = 50.0; //[Hz]
	XL_trans = 0.1; //[pu]
	Z_trans = 0.001; //[pu]
	Im_trans = 0.01; 
	V1_trans = 127.01706; //[kV]
	V2_trans = 35; //[kV]

    //PImodel_2
    double f_2, l_2, rp_2, xp_2, xcp_2, rz_2, xz_2, xcz_2;
    f_2 = 50.0; //[Hz]
    l_2 = 1.0E4; //[m]
    rp_2 = 0.1781598E-4; //[ohn/m]
    xp_2 = 0.31388E-3; //[ohm/m]
    xcp_2 = 273.5448; //[Mohm*m]
    rz_2 = 0.2952201E-3; //[ohm/m]
    xz_2 = 0.1039898E-2; //[ohn/m]
    xcz_2 = 414.1642; //[Mohm*m]

    double *h_Source_pars = { 0 }, *h_Fault_pars = { 0 }, *h_PImodel_1_pars = { 0 }, *h_Trans_sys_pars = { 0 }, *h_PImodel_2_pars = { 0 };
    double *d_Source_pars = { 0 }, *d_Fault_pars = { 0 }, *d_PImodel_1_pars = { 0 }, *d_Trans_sys_pars = { 0 }, *d_PImodel_2_pars = { 0 };

    h_Source_pars = (double*)malloc(4 * sizeof(double));
    h_Fault_pars = (double*)malloc(2 * sizeof(double));
    h_PImodel_1_pars = (double*)malloc(8 * sizeof(double));
    h_Trans_sys_pars = (double*)malloc(7 * sizeof(double));
    h_PImodel_2_pars = (double*)malloc(8 * sizeof(double));

    cudaMalloc((void**)&d_Source_pars, 4 * sizeof(double));
    cudaMalloc((void**)&d_Fault_pars, 2 * sizeof(double));
    cudaMalloc((void**)&d_PImodel_1_pars, 8 * sizeof(double));
    cudaMalloc((void**)&d_Trans_sys_pars, 7 * sizeof(double));
    cudaMalloc((void**)&d_PImodel_2_pars, 8 * sizeof(double));

    h_Source_pars[0] = Vm;
    h_Source_pars[1] = freq;
    h_Source_pars[2] = Rin;
    h_Source_pars[3] = L;

    h_Fault_pars[0] = Ron;
    h_Fault_pars[1] = Roff;

    h_PImodel_1_pars[0] = f_1;
    h_PImodel_1_pars[1] = l_1;
    h_PImodel_1_pars[2] = rp_1;
    h_PImodel_1_pars[3] = xp_1;
    h_PImodel_1_pars[4] = xcp_1;
    h_PImodel_1_pars[5] = rz_1;
    h_PImodel_1_pars[6] = xz_1;
    h_PImodel_1_pars[7] = xcz_1;

    h_Trans_sys_pars[0] = Sb_trans;
    h_Trans_sys_pars[1] = fb_trans;
    h_Trans_sys_pars[2] = XL_trans;
    h_Trans_sys_pars[3] = Z_trans;
    h_Trans_sys_pars[4] = Im_trans;
    h_Trans_sys_pars[5] = V1_trans;
    h_Trans_sys_pars[6] = V2_trans;

    h_PImodel_2_pars[0] = f_2;
    h_PImodel_2_pars[1] = l_2;
    h_PImodel_2_pars[2] = rp_2;
    h_PImodel_2_pars[3] = xp_2;
    h_PImodel_2_pars[4] = xcp_2;
    h_PImodel_2_pars[5] = rz_2;
    h_PImodel_2_pars[6] = xz_2;
    h_PImodel_2_pars[7] = xcz_2;

    cudaMemcpy(d_Source_pars, h_Source_pars, 4 * sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Fault_pars, h_Fault_pars, 2 * sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_PImodel_1_pars, h_PImodel_1_pars, 8 * sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Trans_sys_pars, h_Trans_sys_pars, 7 * sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_PImodel_2_pars, h_PImodel_2_pars, 8 * sizeof(double), cudaMemcpyHostToDevice);

    //////////////////////////////////System Initialization/////////////////////////////////////
    //Source
    Inductor *L_source;
    cudaMalloc((void**)&L_source, 3 * sizeof(Inductor));

    //Fault
    Fault *Fault_1, *Fault_2;
	cudaMalloc((void**)&Fault_1, sizeof(Fault));
    cudaMalloc((void**)&Fault_2, 3 * sizeof(Fault));

    //Transformer
    cudaMalloc((void**)&d_G_trans_sys, 2 * 2 * sizeof(double));
    cudaMalloc((void**)&d_parB_trans_sys, 2 * 2 * sizeof(double));
    cudaMalloc((void**)&d_J_trans_A_sys_d, 2 * 1 * sizeof(double));
    cudaMalloc((void**)&d_J_trans_B_sys_d, 2 * 1 * sizeof(double));
    cudaMalloc((void**)&d_J_trans_C_sys_d, 2 * 1 * sizeof(double));

    //PImodel
    PImodel *PImodel_1, *PImodel_2;
    cudaMalloc((void**)&PImodel_1, sizeof(PImodel));
    cudaMalloc((void**)&PImodel_2, sizeof(PImodel));

    cudaMalloc((void**)&d_Ih_mutual1_1, 3 * sizeof(double));
    cudaMalloc((void**)&d_Ih_mutual2_1, 3 * sizeof(double));
    cudaMalloc((void**)&d_Ih_self1_1, 3 * sizeof(double));
    cudaMalloc((void**)&d_Ih_self2_1, 3 * sizeof(double));
    cudaMalloc((void**)&d_Ih_ABC_1, 3 * sizeof(double));
    
    cudaMalloc((void**)&d_Ih_mutual1_2, 4 * 3 * sizeof(double));
    cudaMalloc((void**)&d_Ih_mutual2_2, 4 * 3 * sizeof(double));
    cudaMalloc((void**)&d_Ih_self1_2, 4 * 3 * sizeof(double));
    cudaMalloc((void**)&d_Ih_self2_2, 4 * 3 * sizeof(double));
    cudaMalloc((void**)&d_Ih_ABC_2, 4 * 3 * sizeof(double));

    //Y Calculation
	//Read A
	//Open file
	FILE* fp = fopen("A.txt", "r");
	
	if(fp == NULL)
	{
		std::cerr << "Cannot open A.txt\n";
        return 1;
	}

	int rows, cols;

	fscanf(fp, "%d %d", &rows, &cols);

	for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            double value;
            fscanf(fp, "%lf", &value);
            A_CAL(i, j) = value;
        }
    }

	fclose(fp);

    cout << A_CAL << endl;

    cudaMalloc((void**)&d_A_CAL_d, sizeof(double) * NODE_NUM * BRANCH_NUM);
    cudaMemcpy(d_A_CAL_d, A_CAL.data(), sizeof(double) * NODE_NUM * BRANCH_NUM, cudaMemcpyHostToDevice);

    cudaMalloc((void**)&d_Yg_CAL_d, sizeof(double) * BRANCH_NUM * BRANCH_NUM);
    cudaMalloc((void**)&d_AYg_CAL_d, sizeof(double) * NODE_NUM * BRANCH_NUM);
    cudaMalloc((void**)&d_Y_CAL_d, sizeof(double) * NODE_NUM * NODE_NUM);
    cudaMalloc((void**)&d_J_branch_CAL_d, sizeof(double) * BRANCH_NUM * 1);
    cudaMalloc((void**)&d_J_CAL_d, sizeof(double) * NODE_NUM * 1);

    cudaMalloc((void**)&d_Y_total_d, sizeof(double) * NODE_NUM * NODE_NUM);
    cudaMalloc((void**)&d_J_total_d, sizeof(double) * NODE_NUM * 1);
    cudaMalloc((void**)&d_Y_total_inv_d, sizeof(double) * NODE_NUM * NODE_NUM);
    cudaMalloc((void**)&d_V_total_d, sizeof(double) * NODE_NUM * 1);

    //P Q
    double *d_Pstation, *d_Qstation, *d_Vstation, *d_Istation;

    cudaMalloc((void**)&d_Pstation, sizeof(double));
    cudaMalloc((void**)&d_Qstation, sizeof(double));
    cudaMalloc((void**)&d_Vstation, sizeof(double) * 3);
    cudaMalloc((void**)&d_Istation, sizeof(double) * 3);

    double *output_his_LPF_Pstation = {0}, double *input_his_LPF_Pstation = {0};
    double *output_his_LPF_Qstation = {0}, double* input_his_LPF_Qstation = {0};

    cudaMalloc((void**)&output_his_LPF_Pstation, sizeof(double));
    cudaMalloc((void**)&input_his_LPF_Pstation, sizeof(double));
    cudaMalloc((void**)&output_his_LPF_Qstation, sizeof(double));
    cudaMalloc((void**)&input_his_LPF_Qstation, sizeof(double));

    ////////////////////////PV Unit Paramter//////////////////////////////
    //DC side
    //PV array
    double Voc_PV, Isc_PV, Vm_PV, Im_PV;
    double Npar_PV, Nser_PV;
	double Tref, Sref;
	double INSO_ini, TEMP_ini, Vpv_ini;
	double par_PV_a, par_PV_b, par_PV_c;

	Voc_PV = 21.7; //[V]
	Isc_PV = 3.35; //[A]
	Vm_PV = 17.4; //[V]
	Im_PV = 3.05; //[A]
	Npar_PV = 656;
	Nser_PV = 29;
	Tref = 25;
	Sref = 1000;
	INSO_ini = 1000; //[W/m^2]
	TEMP_ini = 25; //[Degree] 
	Vpv_ini = 0; //[kV]
	par_PV_a = 0.0025;
	par_PV_b = 0.5;
	par_PV_c = 0.00288;
    
    //PVbreaker and CPVarray
    double Roff_PV = 1e6; //[ohms]
	double Ron_PV = 0.001; //[ohms]
    double CPV = 5000; //[uF]

    //Boost Parameter
    double LBoost = 250e-6; //[H]
    double CBoost = 32000; //[uF]
    double Ron_IGBT = 0.001; //[ohm]
	double Roff_IGBT = 1e8; //[ohm]
    double Ron_Diode = 0.001; //[ohm]
	double Roff_Diode = 1e8; //[ohm]

    //Precharge and Chopper Branches
    double Roff_Chop = 1e8; //[ohm]
	double Ron_Chop = 0.2; //[ohm]
    double Roff_Char = 1e8; //[ohm]
	double Ron_Char = 0.001; //[ohm]

    //VSC Parameter
    double Ron_VSC = 1e-4; //[ohm]
	double Roff_VSC = 1e8; //[ohm]
	double Rsnubber = 0.001; //[ohm]
    double CBus = 32000; //[uF]

    //AC side
    //LC filt Parameter
    double L_LC = 160e-6; //[H]
	double R1_LC = 1e-6; //[ohm]
	double C_LC = 500; //[uF]
	double R2_LC = 0.2; //[ohm]

    //PVunit Breakers
    double Roff_PVunit = 1e4; //[ohms]
	double Ron_PVunit = 1e-3; //[ohms]

    //Transformer
	double Sb_trans_PVunit = 0.333; //[MVA]
	double fb_trans_PVunit = 50; //[Hz]
	double XL_trans_PVunit = 0.15; //[pu]
	double Z_trans_PVunit = 0.001; //[pu]
	double Im_trans_PVunit = 0.004; 
	double V1_trans_PVunit = 0.18187; //[kV]
	double V2_trans_PVunit = 35; //[kV]

    double *h_PVarray_pars = { 0 }, *h_Boost_pars = { 0 }, *h_CPB_pars = { 0 }, *h_VSC_pars = { 0 };
    double *h_LC_pars = { 0 }, *h_PVunit_Bre_pars = { 0 }, *h_Trans_pars = { 0 };
    
    double *d_PVarray_pars = { 0 }, *d_Boost_pars = { 0 }, *d_CPB_pars = { 0 }, *d_VSC_pars = { 0 };
    double *d_LC_pars = { 0 }, *d_PVunit_Bre_pars = { 0 }, *d_Trans_pars = { 0 };

    h_PVarray_pars = (double*)malloc(17 * sizeof(double));
    h_Boost_pars = (double*)malloc(6 * sizeof(double));
    h_CPB_pars = (double*)malloc(4 * sizeof(double));
    h_VSC_pars = (double*)malloc(4 * sizeof(double));
    h_LC_pars = (double*)malloc(4 * sizeof(double));
    h_PVunit_Bre_pars = (double*)malloc(2 * sizeof(double));
    h_Trans_pars = (double*)malloc(7 * sizeof(double));
    
    cudaMalloc((void**)&d_PVarray_pars, 17 * sizeof(double));
    cudaMalloc((void**)&d_Boost_pars, 6 * sizeof(double));
    cudaMalloc((void**)&d_CPB_pars, 4 * sizeof(double));
    cudaMalloc((void**)&d_VSC_pars, 4 * sizeof(double));
    cudaMalloc((void**)&d_LC_pars, 4 * sizeof(double));
    cudaMalloc((void**)&d_PVunit_Bre_pars, 2 * sizeof(double));
    cudaMalloc((void**)&d_Trans_pars, 7 * sizeof(double));

    h_PVarray_pars[0] = Voc_PV;
    h_PVarray_pars[1] = Isc_PV;
	h_PVarray_pars[2] = Vm_PV;
    h_PVarray_pars[3] = Im_PV;
	h_PVarray_pars[4] = Npar_PV;
	h_PVarray_pars[5] = Nser_PV;
	h_PVarray_pars[6] = Tref;
	h_PVarray_pars[7] = Sref;
	h_PVarray_pars[8] = INSO_ini;
	h_PVarray_pars[9] = TEMP_ini; 
	h_PVarray_pars[10] = Vpv_ini;
	h_PVarray_pars[11] = par_PV_a;
	h_PVarray_pars[12] = par_PV_b;
	h_PVarray_pars[13] = par_PV_c;
    h_PVarray_pars[14] = Roff_PV ;
	h_PVarray_pars[15] = Ron_PV;
    h_PVarray_pars[16] = CPV;

    h_Boost_pars[0] = LBoost;
    h_Boost_pars[1] = CBoost;
    h_Boost_pars[2] = Ron_IGBT;
	h_Boost_pars[3] = Roff_IGBT;
    h_Boost_pars[4] = Ron_Diode;
	h_Boost_pars[5] = Roff_Diode;

    h_CPB_pars[0] = Roff_Chop;
	h_CPB_pars[1] = Ron_Chop;
    h_CPB_pars[2] = Roff_Char;
	h_CPB_pars[3] = Ron_Char;

    h_VSC_pars[0] = Ron_VSC;
	h_VSC_pars[1] = Roff_VSC;
	h_VSC_pars[2] = Rsnubber;
    h_VSC_pars[3] = CBus;

    h_LC_pars[0] = L_LC;
	h_LC_pars[1] = R1_LC;
	h_LC_pars[2] = C_LC;
	h_LC_pars[3] = R2_LC;

    h_PVunit_Bre_pars[0] = Roff_PVunit;
    h_PVunit_Bre_pars[1] = Ron_PVunit;

    h_Trans_pars[0] = Sb_trans_PVunit;
	h_Trans_pars[1] = fb_trans_PVunit;
	h_Trans_pars[2] = XL_trans_PVunit;
	h_Trans_pars[3] = Z_trans_PVunit;
	h_Trans_pars[4] = Im_trans_PVunit;
	h_Trans_pars[5] = V1_trans_PVunit;
	h_Trans_pars[6] = V2_trans_PVunit;

    cudaMemcpy(d_PVarray_pars, h_PVarray_pars, 17 * sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Boost_pars, h_Boost_pars, 6 * sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_CPB_pars, h_CPB_pars, 4 * sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_VSC_pars, h_VSC_pars, 4 * sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_LC_pars, h_LC_pars, 4 * sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_PVunit_Bre_pars, h_PVunit_Bre_pars, 2 * sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Trans_pars, h_Trans_pars, 7 * sizeof(double), cudaMemcpyHostToDevice);

    ////////////////////////PV Unit Initialization//////////////////////////////
    //将单个光伏单元中的开关、电容、电感写入一个类
    Switch *Bre_PVunit, *Bre_PV, *Bre_Chop, *Bre_Char, *Boost_IGBT, *Boost_Diode;
    VSC* VSC_average;
    Capacitor* C_PV, * C_Boost, * C_VSC, * C_filt;
    Inductor* L_Boost, * L_filt;

    cudaMalloc((void**)&Bre_PVunit, DIM * 4 * sizeof(Switch));
    cudaMalloc((void**)&Bre_PV, DIM * 4 * sizeof(Switch));
    cudaMalloc((void**)&Bre_Chop, DIM * 4 * sizeof(Switch));
    cudaMalloc((void**)&Bre_Char, DIM * 4 * sizeof(Switch));
    cudaMalloc((void**)&Boost_IGBT, DIM * 4 * sizeof(Switch));
    cudaMalloc((void**)&Boost_Diode, DIM * 4 * sizeof(Switch));

    cudaMalloc((void**)&VSC_average, DIM * 4 * sizeof(VSC));

    cudaMalloc((void**)&C_PV, DIM * 4 * sizeof(Capacitor));
    cudaMalloc((void**)&C_Boost, DIM * 4 * sizeof(Capacitor));
    cudaMalloc((void**)&C_VSC, DIM * 4 * sizeof(Capacitor));
    cudaMalloc((void**)&C_filt, DIM * 4 * 3 * sizeof(Capacitor));

    cudaMalloc((void**)&L_Boost, DIM * 4 * sizeof(Inductor));
    cudaMalloc((void**)&L_filt, DIM * 4 * 3 * sizeof(Inductor));

    //PV arrays
    double *V_PV_V_his = { 0 }, *V_PV_V = { 0 }, *I_PV_A = { 0 }, *V_PV_kV = { 0 }, *I_PV_kA = { 0 };
    double *par_PV_c1 = { 0 }, *par_PV_b1 = { 0 };
    double *Nser = { 0 }, *Npar = { 0 }, * Isc1_PV = { 0 }, * Voc1_PV = { 0 };
    
    cudaMalloc((void**)&I_PV_kA, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&V_PV_kV, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&V_PV_V_his, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&V_PV_V, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&I_PV_A, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&par_PV_c1, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&par_PV_b1, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&Nser, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&Npar, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&Isc1_PV, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&Voc1_PV, sizeof(double) * DIM * 4);

    //Boost
    double* STATUS_his = { 0 };
    cudaMalloc((void**)&STATUS_his, sizeof(double) * 2 * (DIM * 4));

    //LC Filter
    double *G_filt_R1 = { 0 }, *G_filt_R2 = { 0 };
    cudaMalloc((void**)&G_filt_R1, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&G_filt_R2, sizeof(double) * DIM * 4);

    //Transformer
    cudaMalloc((void**)&d_G_trans, (DIM * 4) * (2 * 2) * sizeof(double));
    cudaMalloc((void**)&d_parB_trans, (DIM * 4) * (2 * 2) * sizeof(double));
    cudaMalloc((void**)&d_J_trans_A_PVunit_d, (DIM * 4) * 2 * sizeof(double));
    cudaMalloc((void**)&d_J_trans_B_PVunit_d, (DIM * 4) * 2 * sizeof(double));
    cudaMalloc((void**)&d_J_trans_C_PVunit_d, (DIM * 4) * 2 * sizeof(double));

    cudaMalloc((void**)&I_port_A_PVunit, (DIM * 4) * 2 * sizeof(double));
    cudaMalloc((void**)&I_port_B_PVunit, (DIM * 4) * 2 * sizeof(double));
    cudaMalloc((void**)&I_port_C_PVunit, (DIM * 4) * 2 * sizeof(double));

    ///////////////////////////Matrix Parameter////////////////////////////////////
    //读入数据
    FILE* fp1;
    FILE* fp2;
    errno_t err1 = fopen_s(&fp1, "A1.txt", "r");
    if (err1 != 0) {
        printf("无法打开文件1\n");
        return -1;
    }
    errno_t err2 = fopen_s(&fp2, "A2.txt", "r");
    if (err2 != 0) {
        printf("无法打开文件2\n");
        return -1;
    }

    int rows1, cols1, rows2, cols2;
    fscanf(fp1, "%d %d", &rows1, &cols1);

    for (int i = 0; i < rows1; i++) {
        for (int j = 0; j < cols1; j++) {
            double value;
            fscanf(fp1, "%lf", &value);
            A1(i, j) = value;
        }
    }

    fscanf(fp2, "%d %d", &rows2, &cols2);

    for (int i = 0; i < rows2; i++) {
        for (int j = 0; j < cols2; j++) {
            double value;
            fscanf(fp2, "%lf", &value);
            A2(i, j) = value;
        }
    }

    fclose(fp1);
    fclose(fp2);

    //Y1 calculation (AC side)
    double **d_A1, **d_Yg1, **d_AYg1, **d_Y1, **d_J_branch1, **d_J1;
    double **d_Y11, **d_Y12, **d_Y21, **d_Y22, **d_Y22_inv, **d_YY, **d_YYY, **d_Yeq;
    double **d_Jex, **d_Jin, **d_YYJ, **d_Jeq;
    double **d_Vin, **d_Vex, **d_YV, **d_JYV;

    //A1
    d_A1 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_A1[i], sizeof(double) * NODE_NUM_AC * BRANCH_NUM_AC);
        cudaMemcpy(d_A1[i], A1.data(), sizeof(double) * NODE_NUM_AC * BRANCH_NUM_AC, cudaMemcpyHostToDevice);
    }

    cudaMalloc(&d_A1_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_A1_d, d_A1, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Yg1
    d_Yg1 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Yg1[i], sizeof(double) * BRANCH_NUM_AC * BRANCH_NUM_AC);
    }

    cudaMalloc(&d_Yg1_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Yg1_d, d_Yg1, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //AYg1
    d_AYg1 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_AYg1[i], sizeof(double) * NODE_NUM_AC * BRANCH_NUM_AC);
    }

    cudaMalloc(&d_AYg1_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_AYg1_d, d_AYg1, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Y1
    d_Y1 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Y1[i], sizeof(double) * NODE_NUM_AC * NODE_NUM_AC);
    }

    cudaMalloc(&d_Y1_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Y1_d, d_Y1, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //J_branch1
    d_J_branch1 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_J_branch1[i], sizeof(double) * BRANCH_NUM_AC * 1);
    }

    cudaMalloc(&d_J_branch1_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_J_branch1_d, d_J_branch1, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //J1
    d_J1 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_J1[i], sizeof(double) * NODE_NUM_AC * 1);
    }

    cudaMalloc(&d_J1_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_J1_d, d_J1, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //节点消去部分(AC side)
    //Y11
    d_Y11 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Y11[i], sizeof(double) * 3 * 3);
    }

    cudaMalloc(&d_Y11_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Y11_d, d_Y11, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Y12
    d_Y12 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Y12[i], sizeof(double) * 3 * (NODE_NUM_AC - 3));
    }

    cudaMalloc(&d_Y12_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Y12_d, d_Y12, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Y21
    d_Y21 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Y21[i], sizeof(double) * (NODE_NUM_AC - 3) * 3);
    }

    cudaMalloc(&d_Y21_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Y21_d, d_Y21, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Y22
    d_Y22 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Y22[i], sizeof(double) * (NODE_NUM_AC - 3) * (NODE_NUM_AC - 3));
    }

    cudaMalloc(&d_Y22_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Y22_d, d_Y22, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);
    
    //Y22_inv
    d_Y22_inv = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Y22_inv[i], sizeof(double) * (NODE_NUM_AC - 3) * (NODE_NUM_AC - 3));
    }

    cudaMalloc(&d_Y22_inv_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Y22_inv_d, d_Y22_inv, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //YY
    d_YY = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_YY[i], sizeof(double) * 3 * (NODE_NUM_AC - 3));
    }

    cudaMalloc(&d_YY_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_YY_d, d_YY, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //YYY
    d_YYY = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_YYY[i], sizeof(double) * 3 * 3);
    }

    cudaMalloc(&d_YYY_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_YYY_d, d_YYY, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Yeq
    d_Yeq = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Yeq[i], sizeof(double) * 3 * 3);
    }

    cudaMalloc(&d_Yeq_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Yeq_d, d_Yeq, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Jex
    d_Jex = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Jex[i], sizeof(double) * 3 * 1);
    }

    cudaMalloc(&d_Jex_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Jex_d, d_Jex, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Jin
    d_Jin = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Jin[i], sizeof(double) * (NODE_NUM_AC - 3) * 1);
    }

    cudaMalloc(&d_Jin_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Jin_d, d_Jin, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //YYJ
    d_YYJ = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_YYJ[i], sizeof(double) * 3 * 1);
    }

    cudaMalloc(&d_YYJ_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_YYJ_d, d_YYJ, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Jeq
    d_Jeq = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Jeq[i], sizeof(double) * 3 * 1);
    }

    cudaMalloc(&d_Jeq_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Jeq_d, d_Jeq, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Vin
    d_Vin = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Vin[i], sizeof(double) * (NODE_NUM_AC - 3) * 1);
    }

    cudaMalloc(&d_Vin_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Vin_d, d_Vin, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Vex
    d_Vex = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Vex[i], sizeof(double) * 3 * 1);
    }

    cudaMalloc(&d_Vex_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Vex_d, d_Vex, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //YV
    d_YV = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_YV[i], sizeof(double) * (NODE_NUM_AC - 3) * 1);
    }

    cudaMalloc(&d_YV_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_YV_d, d_YV, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //JYV
    d_JYV = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_JYV[i], sizeof(double) * (NODE_NUM_AC - 3) * 1);
    }

    cudaMalloc(&d_JYV_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_JYV_d, d_JYV, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Y2 calculation (DC side)
    double **d_A2, **d_Yg2, **d_AYg2, **d_Y2;
    double **d_J_branch2, **d_J2;
    double **d_V_DC, **d_Y2_inv;

    //A2
    d_A2 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_A2[i], sizeof(double) * NODE_NUM_DC * BRANCH_NUM_DC);
        cudaMemcpy(d_A2[i], A2.data(), sizeof(double) * NODE_NUM_DC * BRANCH_NUM_DC, cudaMemcpyHostToDevice);
    }

    cudaMalloc(&d_A2_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_A2_d, d_A2, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Yg2
    d_Yg2 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Yg2[i], sizeof(double) * BRANCH_NUM_DC * BRANCH_NUM_DC);
    }

    cudaMalloc(&d_Yg2_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Yg2_d, d_Yg2, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //AYg2
    d_AYg2 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_AYg2[i], sizeof(double) * NODE_NUM_DC * BRANCH_NUM_DC);
    }

    cudaMalloc(&d_AYg2_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_AYg2_d, d_AYg2, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Y2
    d_Y2 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_Y2[i], sizeof(double) * NODE_NUM_DC * NODE_NUM_DC);
    }

    cudaMalloc(&d_Y2_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Y2_d, d_Y2, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //J_branch2
    d_J_branch2 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_J_branch2[i], sizeof(double) * BRANCH_NUM_DC * 1);
    }

    cudaMalloc(&d_J_branch2_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_J_branch2_d, d_J_branch2, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //J2
    d_J2 = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_J2[i], sizeof(double) * NODE_NUM_DC * 1);
    }

    cudaMalloc(&d_J2_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_J2_d, d_J2, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //VDC
    d_V_DC = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; i++)
    {
        cudaMalloc((void**)&d_V_DC[i], sizeof(double) * NODE_NUM_DC * 1);
    }

    cudaMalloc(&d_V_DC_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_V_DC_d, d_V_DC, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    //Y2_inv
    d_Y2_inv = new double* [DIM * 4];
    for (int i = 0; i < DIM * 4; ++i) 
    {
        cudaMalloc(&d_Y2_inv[i], NODE_NUM_DC * NODE_NUM_DC * sizeof(double));
    }

    cudaMalloc(&d_Y2_inv_d, sizeof(double*) * DIM * 4);
    cudaMemcpy(d_Y2_inv_d, d_Y2_inv, sizeof(double*) * DIM * 4, cudaMemcpyHostToDevice);

    ////////////////////////Control Paramter//////////////////////////////
    //initial control
    //initialize the pointer to NULL
    double *h_fsys = { 0 };
    double *d_fsys = { 0 };
    
    h_fsys = (double*)malloc(sizeof(double));
    cudaMalloc((void**)&d_fsys, sizeof(double));

    *h_fsys = 50;
    cudaMemcpy(d_fsys, h_fsys, sizeof(double), cudaMemcpyHostToDevice);

    int *d_Dblk = { 0 }, *d_SW_Char = { 0 };
    
    cudaMalloc((void**)&d_Dblk, sizeof(int) * 4 * DIM);
    cudaMalloc((void**)&d_SW_Char, sizeof(int) * 4 * DIM);

    //start control and Parameters set
    double *h_VdcPURef = { 0 }, *h_QPURef = { 0 }, *h_VacPURef = { 0 };
    int *h_VacOrQ = { 0 }, *h_Har3ctrl = { 0 }, *h_FrzPI = { 0 }, *h_VdcCon = { 0 };

    double *d_VdcPURef = { 0 }, *d_QPURef = { 0 }, *d_VacPURef = { 0 };
    int *d_VacOrQ = { 0 }, *d_Har3ctrl = { 0 }, *d_FrzPI = { 0 }, *d_VdcCon = { 0 };

    h_VdcPURef = (double*)malloc(sizeof(double));
    h_QPURef = (double*)malloc(sizeof(double));
    h_VacPURef = (double*)malloc(sizeof(double));
    h_VacOrQ = (int*)malloc(sizeof(int));
    h_Har3ctrl = (int*)malloc(sizeof(int));
    h_FrzPI = (int*)malloc(sizeof(int));
    h_VdcCon = (int*)malloc(sizeof(int));

    cudaMalloc((void**)&d_VdcPURef, sizeof(double));
    cudaMalloc((void**)&d_QPURef, sizeof(double));
    cudaMalloc((void**)&d_VacPURef, sizeof(double));
    cudaMalloc((void**)&d_VacOrQ, sizeof(int));
    cudaMalloc((void**)&d_Har3ctrl, sizeof(int));
    cudaMalloc((void**)&d_FrzPI, sizeof(int));
    cudaMalloc((void**)&d_VdcCon, sizeof(int));

    *h_VdcPURef = 1;
    *h_QPURef = -0.14;
    *h_VacPURef = 1;
    *h_VacOrQ = 1; //开关信号1/0
    *h_Har3ctrl = 1; //开关信号1/0
    *h_FrzPI = 0;
    *h_VdcCon = 0; //开关信号1/0

    cudaMemcpy(d_VdcPURef, h_VdcPURef, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_QPURef, h_QPURef, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_VacPURef, h_VacPURef, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_VacOrQ, h_VacOrQ, sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Har3ctrl, h_Har3ctrl, sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_FrzPI, h_FrzPI, sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_VdcCon, h_VdcCon, sizeof(int), cudaMemcpyHostToDevice);

    //Vsc control
    double *h_VLLrated = { 0 }, *h_VdcRated = { 0 }, *h_PoRated = { 0 }, *h_VpvRated = { 0 };
    double *h_Lcon1 = { 0 }, *h_fsw1 = { 0 }, *h_fsw2 = { 0 }, *h_Vtr2Pri = { 0 }, *h_Vtr2Sec = { 0 };

    double *d_VLLrated = { 0 }, *d_VdcRated = { 0 }, *d_PoRated = { 0 }, *d_VpvRated = { 0 };
    double *d_Lcon1 = { 0 }, *d_fsw1 = { 0 }, *d_fsw2 = { 0 }, *d_Vtr2Pri = { 0 }, *d_Vtr2Sec = { 0 };

    h_VLLrated = (double*)malloc(sizeof(double));
    h_VdcRated = (double*)malloc(sizeof(double));
    h_PoRated = (double*)malloc(sizeof(double));
    h_VpvRated = (double*)malloc(sizeof(double));
    h_Lcon1 = (double*)malloc(sizeof(double));
    h_fsw1 = (double*)malloc(sizeof(double));
    h_fsw2 = (double*)malloc(sizeof(double));
    h_Vtr2Pri = (double*)malloc(sizeof(double));
    h_Vtr2Sec = (double*)malloc(sizeof(double));

    cudaMalloc((void**)&d_VLLrated, sizeof(double));
    cudaMalloc((void**)&d_VdcRated, sizeof(double));
    cudaMalloc((void**)&d_PoRated, sizeof(double));
    cudaMalloc((void**)&d_VpvRated, sizeof(double));
    cudaMalloc((void**)&d_Lcon1, sizeof(double));
    cudaMalloc((void**)&d_fsw1, sizeof(double));
    cudaMalloc((void**)&d_fsw2, sizeof(double));
    cudaMalloc((void**)&d_Vtr2Pri, sizeof(double));
    cudaMalloc((void**)&d_Vtr2Sec, sizeof(double));

    *h_VLLrated = 35.0;
    *h_VdcRated = 0.8;
    *h_PoRated = 1.0;
    *h_VpvRated = 0.5;
    *h_Lcon1 = 0.00016;
    *h_fsw1 = 2000.0;
    *h_fsw2 = 2000.0;
    *h_Vtr2Pri = 35.0; //VLL for D windings
    *h_Vtr2Sec = 0.18187; //Vph for yn windings

    cudaMemcpy(d_VLLrated, h_VLLrated, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_VdcRated, h_VdcRated, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_PoRated, h_PoRated, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_VpvRated, h_VpvRated, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Lcon1, h_Lcon1, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_fsw1, h_fsw1, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_fsw2, h_fsw2, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Vtr2Pri, h_Vtr2Pri, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Vtr2Sec, h_Vtr2Sec, sizeof(double), cudaMemcpyHostToDevice);

    //Control PI Parameter
    double *h_fmax1 = { 0 }, *h_IdMax = { 0 }, *h_IdMin = { 0 };
    double *d_fmax1 = { 0 }, *d_IdMax = { 0 }, *d_IdMin = { 0 };

    h_fmax1 = (double*)malloc(sizeof(double));
    h_IdMax = (double*)malloc(sizeof(double));
    h_IdMin = (double*)malloc(sizeof(double));

    cudaMalloc((void**)&d_fmax1, sizeof(double));
    cudaMalloc((void**)&d_IdMax, sizeof(double));
    cudaMalloc((void**)&d_IdMin, sizeof(double));
    
    *h_fmax1 = 55;
    *h_IdMax = 1.5;
    *h_IdMin = -1.5;

    cudaMemcpy(d_fmax1, h_fmax1, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_IdMax, h_IdMax, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_IdMin, h_IdMin, sizeof(double), cudaMemcpyHostToDevice);

    //Control
    //vsc
    double *h_G_VSC = { 0 }, *h_T_VSC = { 0 };
    double *d_G_VSC = { 0 }, *d_T_VSC = { 0 };

    h_G_VSC = (double*)malloc(sizeof(double));
    h_T_VSC = (double*)malloc(sizeof(double));

    cudaMalloc((void**)&d_G_VSC, sizeof(double));
    cudaMalloc((void**)&d_T_VSC, sizeof(double));

    *h_G_VSC = 1.0;
    *h_T_VSC = 0.01;

    cudaMemcpy(d_G_VSC, h_G_VSC, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_T_VSC, h_T_VSC, sizeof(double), cudaMemcpyHostToDevice);
    
    //Boost Control
    double *h_Strtup = { 0 }, *h_Samp = { 0 }, *h_Vstep = { 0 }, *h_Sd = { 0 };
    double *h_G1_Boost = { 0 }, *h_T1_Boost = { 0 }, *h_G2_Boost = { 0 }, *h_T2_Boost = { 0 }, *h_G3_Boost = { 0 }, *h_T3_Boost = { 0 };
    int *h_signal = { 0 };
    double *h_Max_Boost = { 0 }, *h_Min_Boost = { 0 }, *h_Kp_Boost = { 0 }, *h_Ti_Boost = { 0 };

    double *d_Strtup = { 0 }, *d_Samp = { 0 }, *d_Vstep = { 0 }, *d_Sd = { 0 };
    double *d_G1_Boost = { 0 }, *d_T1_Boost = { 0 }, *d_G2_Boost = { 0 }, *d_T2_Boost = { 0 }, *d_G3_Boost = { 0 }, *d_T3_Boost = { 0 };
    int *d_signal = { 0 };
    double *d_Max_Boost = { 0 }, *d_Min_Boost = { 0 }, *d_Kp_Boost = { 0 }, *d_Ti_Boost = { 0 };

    h_Strtup = (double*)malloc(sizeof(double));
    h_Samp = (double*)malloc(sizeof(double));
    h_Vstep = (double*)malloc(sizeof(double));
    h_Sd = (double*)malloc(sizeof(double));
    h_G1_Boost = (double*)malloc(sizeof(double));
    h_T1_Boost = (double*)malloc(sizeof(double));
    h_G2_Boost = (double*)malloc(sizeof(double));
    h_T2_Boost = (double*)malloc(sizeof(double));
    h_G3_Boost = (double*)malloc(sizeof(double));
    h_T3_Boost = (double*)malloc(sizeof(double));
    h_Max_Boost = (double*)malloc(sizeof(double));
    h_Min_Boost = (double*)malloc(sizeof(double));
    h_Kp_Boost = (double*)malloc(sizeof(double));
    h_Ti_Boost = (double*)malloc(sizeof(double));

    h_signal = (int*)malloc(sizeof(int));
    
    cudaMalloc((void**)&d_Strtup, sizeof(double));
    cudaMalloc((void**)&d_Samp, sizeof(double));
    cudaMalloc((void**)&d_Vstep, sizeof(double));
    cudaMalloc((void**)&d_Sd, sizeof(double));
    cudaMalloc((void**)&d_G1_Boost, sizeof(double));
    cudaMalloc((void**)&d_T1_Boost, sizeof(double));
    cudaMalloc((void**)&d_G2_Boost, sizeof(double));
    cudaMalloc((void**)&d_T2_Boost, sizeof(double));
    cudaMalloc((void**)&d_G3_Boost, sizeof(double));
    cudaMalloc((void**)&d_T3_Boost, sizeof(double));
    cudaMalloc((void**)&d_Max_Boost, sizeof(double));
    cudaMalloc((void**)&d_Min_Boost, sizeof(double));
    cudaMalloc((void**)&d_Kp_Boost, sizeof(double));
    cudaMalloc((void**)&d_Ti_Boost, sizeof(double));

    cudaMalloc((void**)&d_signal, sizeof(int));

    *h_Strtup = 2.5;
    *h_Samp = 50;
    *h_Vstep = 10;
    *h_Sd = 0.02;
    *h_G1_Boost = 1.0;
    *h_T1_Boost = 0.05;
    *h_G2_Boost = 1.0;
    *h_T2_Boost = 0.05;
    *h_G3_Boost = 1.0;
    *h_T3_Boost = 0.01;
    *h_signal = 1;
    *h_Max_Boost = 0.95;
    *h_Min_Boost = 0;
    *h_Kp_Boost = 0.5;
    *h_Ti_Boost = 0.05;

    cudaMemcpy(d_Strtup, h_Strtup, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Samp, h_Samp, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Vstep, h_Vstep, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Sd, h_Sd, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_G1_Boost, h_G1_Boost, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_T1_Boost, h_T1_Boost, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_G2_Boost, h_G2_Boost, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_T2_Boost, h_T2_Boost, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_G3_Boost, h_G3_Boost, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_T3_Boost, h_T3_Boost, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Max_Boost, h_Max_Boost, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Min_Boost, h_Min_Boost, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Kp_Boost, h_Kp_Boost, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Ti_Boost, h_Ti_Boost, sizeof(double), cudaMemcpyHostToDevice);

    cudaMemcpy(d_signal, h_signal, sizeof(int), cudaMemcpyHostToDevice);

    //Voltage Loop
    double *h_Kp = { 0 }, *h_Kp2 = { 0 }, *h_ULIM = { 0 }, *h_LLIM = { 0 };
    double *h_MAXOUT = { 0 }, *h_MINOUT = { 0 }, *h_MAXOUT2 = { 0 }, *h_MINOUT2 = { 0 };
    double *h_PI_T1 = { 0 }, *h_PI_T2 = { 0 };

    double *d_Kp = { 0 }, *d_Kp2 = { 0 }, *d_ULIM = { 0 }, *d_LLIM = { 0 };
    double *d_MAXOUT = { 0 }, *d_MINOUT = { 0 }, *d_MAXOUT2 = { 0 }, *d_MINOUT2 = { 0 };
    double *d_PI_T1 = { 0 }, *d_PI_T2 = { 0 };

    h_Kp = (double*)malloc(sizeof(double));
    h_Kp2 = (double*)malloc(sizeof(double));
    h_ULIM = (double*)malloc(sizeof(double));
    h_LLIM = (double*)malloc(sizeof(double));
    h_MAXOUT = (double*)malloc(sizeof(double));
    h_MINOUT = (double*)malloc(sizeof(double));
    h_MAXOUT2 = (double*)malloc(sizeof(double));
    h_MINOUT2 = (double*)malloc(sizeof(double));
    h_PI_T1 = (double*)malloc(sizeof(double));
    h_PI_T2 = (double*)malloc(sizeof(double));
    
    cudaMalloc((void**)&d_Kp, sizeof(double));
    cudaMalloc((void**)&d_Kp2, sizeof(double));
    cudaMalloc((void**)&d_ULIM, sizeof(double));
    cudaMalloc((void**)&d_LLIM, sizeof(double));
    cudaMalloc((void**)&d_MAXOUT, sizeof(double));
    cudaMalloc((void**)&d_MINOUT, sizeof(double));
    cudaMalloc((void**)&d_MAXOUT2, sizeof(double));
    cudaMalloc((void**)&d_MINOUT2, sizeof(double));
    cudaMalloc((void**)&d_PI_T1, sizeof(double));
    cudaMalloc((void**)&d_PI_T2, sizeof(double));

    *h_Kp = 1;
    *h_Kp2 = 1;
    *h_ULIM = 10;
    *h_LLIM = 10;
    *h_MAXOUT = 1.5;
    *h_MINOUT = -1.5;
    *h_MAXOUT2 = 1;
    *h_MINOUT2 = -1;
    *h_PI_T1 = 0.02;
    *h_PI_T2 = 0.02;

    cudaMemcpy(d_Kp, h_Kp, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Kp2, h_Kp2, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_ULIM, h_ULIM, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_LLIM, h_LLIM, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_MAXOUT, h_MAXOUT, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_MINOUT, h_MINOUT, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_MAXOUT2, h_MAXOUT2, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_MINOUT2, h_MINOUT2, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_PI_T1, h_PI_T1, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_PI_T2, h_PI_T2, sizeof(double), cudaMemcpyHostToDevice);

    //Current Loop
    double *h_Mmax_current = { 0 }, *h_Mmin_current = { 0 }, *h_MaxLimit1_current = { 0 }, *h_MinLimit1_current = { 0 }, *h_MaxLimit2_current = { 0 }, *h_MinLimit2_current = { 0 };
    double *h_Kp1_current = { 0 }, *h_Kp2_current = { 0 }, *h_Ti1_current = { 0 }, *h_Ti2_current = { 0 };
    
    double *d_Mmax_current = { 0 }, *d_Mmin_current = { 0 }, *d_MaxLimit1_current = { 0 }, *d_MinLimit1_current = { 0 }, *d_MaxLimit2_current = { 0 }, *d_MinLimit2_current = { 0 };
    double *d_Kp1_current = { 0 }, *d_Kp2_current = { 0 }, *d_Ti1_current = { 0 }, *d_Ti2_current = { 0 };

    h_Mmax_current = (double*)malloc(sizeof(double));
    h_Mmin_current = (double*)malloc(sizeof(double));
    h_MaxLimit1_current = (double*)malloc(sizeof(double));
    h_MinLimit1_current = (double*)malloc(sizeof(double));
    h_MaxLimit2_current = (double*)malloc(sizeof(double));
    h_MinLimit2_current = (double*)malloc(sizeof(double));
    h_Kp1_current = (double*)malloc(sizeof(double));
    h_Kp2_current = (double*)malloc(sizeof(double));
    h_Ti1_current = (double*)malloc(sizeof(double));
    h_Ti2_current = (double*)malloc(sizeof(double));
    
    cudaMalloc((void**)&d_Mmax_current, sizeof(double));
    cudaMalloc((void**)&d_Mmin_current, sizeof(double));
    cudaMalloc((void**)&d_MaxLimit1_current, sizeof(double));
    cudaMalloc((void**)&d_MinLimit1_current, sizeof(double));
    cudaMalloc((void**)&d_MaxLimit2_current, sizeof(double));
    cudaMalloc((void**)&d_MinLimit2_current, sizeof(double));
    cudaMalloc((void**)&d_Kp1_current, sizeof(double));
    cudaMalloc((void**)&d_Kp2_current, sizeof(double));
    cudaMalloc((void**)&d_Ti1_current, sizeof(double));
    cudaMalloc((void**)&d_Ti2_current, sizeof(double));

    *h_Mmax_current = 1.25;
    *h_Mmin_current = 0;
    *h_MaxLimit1_current = 1;
    *h_MinLimit1_current = -1;
    *h_MaxLimit2_current = 1;
    *h_MinLimit2_current = -1;
    *h_Kp1_current = 0.5;
    *h_Kp2_current = 0.5;
    *h_Ti1_current = 0.25;
    *h_Ti2_current = 0.1;

    cudaMemcpy(d_Mmax_current, h_Mmax_current, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Mmin_current, h_Mmin_current, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_MaxLimit1_current, h_MaxLimit1_current, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_MinLimit1_current, h_MinLimit1_current, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_MaxLimit2_current, h_MaxLimit2_current, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_MinLimit2_current, h_MinLimit2_current, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Kp1_current, h_Kp1_current, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Kp2_current, h_Kp2_current, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Ti1_current, h_Ti1_current, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Ti2_current, h_Ti2_current, sizeof(double), cudaMemcpyHostToDevice);

    //PLL & DQ
    double *h_Kp_PLL = { 0 }, *h_Ti_PLL = { 0 }, *h_Base_Vpll = { 0 }, *h_fbase = { 0 };
    double *h_fmaxPU = { 0 }, *h_fmin1 = { 0 }, *h_fminPU = { 0 };
    double *h_Kp_realpole = { 0 }, *h_Ti_realpole_L = { 0 }, *h_Ti_realpole_R = { 0 };
    double *h_theta1 = { 0 }, *h_theta2 = { 0 };
    
    double *d_Kp_PLL = { 0 }, *d_Ti_PLL = { 0 }, *d_Base_Vpll = { 0 }, *d_fbase = { 0 };
    double *d_fmaxPU = { 0 }, *d_fmin1 = { 0 }, *d_fminPU = { 0 };
    double *d_Kp_realpole = { 0 }, *d_Ti_realpole_L = { 0 }, *d_Ti_realpole_R = { 0 };
    double *d_theta1 = { 0 }, *d_theta2 = { 0 };

    h_Kp_PLL = (double*)malloc(sizeof(double));
    h_Ti_PLL = (double*)malloc(sizeof(double));
    h_Base_Vpll = (double*)malloc(sizeof(double));
    h_fbase = (double*)malloc(sizeof(double));
    
    h_fmaxPU = (double*)malloc(sizeof(double));
    h_fmin1 = (double*)malloc(sizeof(double));
    h_fminPU = (double*)malloc(sizeof(double));
    h_Kp_realpole = (double*)malloc(sizeof(double));
    h_Ti_realpole_L = (double*)malloc(sizeof(double));
    h_Ti_realpole_R = (double*)malloc(sizeof(double));
    h_theta1 = (double*)malloc(sizeof(double));
    h_theta2 = (double*)malloc(sizeof(double));

    cudaMalloc((void**)&d_Kp_PLL, sizeof(double));
    cudaMalloc((void**)&d_Ti_PLL, sizeof(double));
    cudaMalloc((void**)&d_Base_Vpll, sizeof(double));
    cudaMalloc((void**)&d_fbase, sizeof(double));

    cudaMalloc((void**)&d_fmaxPU, sizeof(double));
    cudaMalloc((void**)&d_fmin1, sizeof(double));
    cudaMalloc((void**)&d_fminPU, sizeof(double));
    cudaMalloc((void**)&d_Kp_realpole, sizeof(double));
    cudaMalloc((void**)&d_Ti_realpole_L, sizeof(double));
    cudaMalloc((void**)&d_Ti_realpole_R, sizeof(double));
    cudaMalloc((void**)&d_theta1, sizeof(double));
    cudaMalloc((void**)&d_theta2, sizeof(double));

    *h_Kp_PLL = 5;
    *h_Ti_PLL = 0.01;
    *h_Base_Vpll = 1;
    *h_fbase = 50;

    *h_fmaxPU = *h_fmax1 / *h_fsys;
    *h_fmin1 = 45;
    *h_fminPU = *h_fmin1 / *h_fsys;
    *h_Kp_realpole = 1.0;
    *h_Ti_realpole_L = 0.01;
    *h_Ti_realpole_R = 0.001;
    *h_theta1 = 0;
    *h_theta2 = 0;

    cudaMemcpy(d_Kp_PLL, h_Kp_PLL, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Ti_PLL, h_Ti_PLL, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Base_Vpll, h_Base_Vpll, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_fbase, h_fbase, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_fmaxPU, h_fmaxPU, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_fmin1, h_fmin1, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_fminPU, h_fminPU, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Kp_realpole, h_Kp_realpole, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Ti_realpole_L, h_Ti_realpole_L, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Ti_realpole_R, h_Ti_realpole_R, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_theta1, h_theta1, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_theta2, h_theta2, sizeof(double), cudaMemcpyHostToDevice);

    //3rd Harmonic
    int *h_Ctrl = { 0 };
    int *d_Ctrl = { 0 };

    h_Ctrl = (int*)malloc(sizeof(int));

    cudaMalloc((void**)&d_Ctrl, sizeof(int));

    *h_Ctrl = 0;

    cudaMemcpy(d_Ctrl, h_Ctrl, sizeof(int), cudaMemcpyHostToDevice);
    
    //PWM
    double *h_max1_pwm = { 0 }, *h_min1_pwm = { 0 }, *h_max2_pwm = { 0 }, *h_min2_pwm = { 0 };
    double *d_max1_pwm = { 0 }, *d_min1_pwm = { 0 }, *d_max2_pwm = { 0 }, *d_min2_pwm = { 0 };

    h_max1_pwm = (double*)malloc(sizeof(double));
    h_min1_pwm = (double*)malloc(sizeof(double));
    h_max2_pwm = (double*)malloc(sizeof(double));
    h_min2_pwm = (double*)malloc(sizeof(double));

    cudaMalloc((void**)&d_max1_pwm, sizeof(double));
    cudaMalloc((void**)&d_min1_pwm, sizeof(double));
    cudaMalloc((void**)&d_max2_pwm, sizeof(double));
    cudaMalloc((void**)&d_min2_pwm, sizeof(double));

    *h_max1_pwm = 1;
    *h_min1_pwm = -1;
    *h_max2_pwm = 1;
    *h_min2_pwm = 0;

    cudaMemcpy(d_max1_pwm, h_max1_pwm, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_min1_pwm, h_min1_pwm, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_max2_pwm, h_max2_pwm, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_min2_pwm, h_min2_pwm, sizeof(double), cudaMemcpyHostToDevice);
    
    //Start
    double *h_G_start = { 0 }, *h_T_start = { 0 }, *h_comparator1_H = { 0 }, *h_comparator1_L = { 0 }, *h_comparator2_H = { 0 }, *h_comparator2_L = { 0 };
    double *d_G_start = { 0 }, *d_T_start = { 0 }, *d_comparator1_H = { 0 }, *d_comparator1_L = { 0 }, *d_comparator2_H = { 0 }, *d_comparator2_L = { 0 };

    h_G_start = (double*)malloc(sizeof(double));
    h_T_start = (double*)malloc(sizeof(double));
    h_comparator1_H = (double*)malloc(sizeof(double));
    h_comparator1_L = (double*)malloc(sizeof(double));
    h_comparator2_H = (double*)malloc(sizeof(double));
    h_comparator2_L = (double*)malloc(sizeof(double));
    
    cudaMalloc((void**)&d_G_start, sizeof(double));
    cudaMalloc((void**)&d_T_start, sizeof(double));
    cudaMalloc((void**)&d_comparator1_H, sizeof(double));
    cudaMalloc((void**)&d_comparator1_L, sizeof(double));
    cudaMalloc((void**)&d_comparator2_H, sizeof(double));
    cudaMalloc((void**)&d_comparator2_L, sizeof(double));

    *h_G_start = 1.0;
    *h_T_start = 0.001;
    *h_comparator1_H = 1;
    *h_comparator1_L = 0;
    *h_comparator2_H = 1;
    *h_comparator2_L = 0;

    cudaMemcpy(d_G_start, h_G_start, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_T_start, h_T_start, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_comparator1_H, h_comparator1_H, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_comparator1_L, h_comparator1_L, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_comparator2_H, h_comparator2_H, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_comparator2_L, h_comparator2_L, sizeof(double), cudaMemcpyHostToDevice);

    //Defination
    double *h_comparator_1_B = { 0 }, *h_comparator_2_A = { 0 }, *h_signal_A = { 0 }, *h_Park0 = { 0 };
    double *d_comparator_1_B = { 0 }, *d_comparator_2_A = { 0 }, *d_signal_A = { 0 }, *d_Park0 = { 0 };

    h_comparator_1_B = (double*)malloc(sizeof(double));
    h_comparator_2_A = (double*)malloc(sizeof(double));
    h_signal_A = (double*)malloc(sizeof(double));
    h_Park0 = (double*)malloc(sizeof(double));

    cudaMalloc((void**)&d_comparator_1_B, sizeof(double));
    cudaMalloc((void**)&d_comparator_2_A, sizeof(double));
    cudaMalloc((void**)&d_signal_A, sizeof(double));
    cudaMalloc((void**)&d_Park0, sizeof(double));

    *h_comparator_1_B = 1.4;
    *h_comparator_2_A = 0.98;
    *h_signal_A = 0.51;
    *h_Park0 = 0.0;

    cudaMemcpy(d_comparator_1_B, h_comparator_1_B, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_comparator_2_A, h_comparator_2_A, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_signal_A, h_signal_A, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_Park0, h_Park0, sizeof(double), cudaMemcpyHostToDevice);

    //Initialization
    double *VdcBusPU= { 0 }, *VdcBus = {0};
    double *Vcon1A = { 0 }, *Vcon1B = { 0 }, *Vcon1C = { 0 };
    double *Icon1A = { 0 }, *Icon1B = { 0 }, *Icon1C = { 0 }, *Icon1A1 = { 0 }, *Icon1B1 = { 0 }, *Icon1C1 = { 0 };
    double *Pcon1filt = { 0 }, *Qcon1filt = { 0 }, *VacRMSfilt = { 0 }, *Ipv = { 0 }, *VpvBus = { 0 };
    double * VdcPreCh = { 0 };
    int *MPPT_ctrl = { 0 }, *Ta_1 = { 0 }, *Ta_4 = { 0 }, *Tb_3 = { 0 }, *Tb_6 = { 0 }, *Tc_5 = { 0 }, *Tc_2 = { 0 };
    double *MPPT_ctrl_timetag = { 0 }, *Ta_1_timetag = { 0 }, *Ta_4_timetag = { 0 }, *Tb_3_timetag = { 0 }, *Tb_6_timetag = { 0 }, *Tc_5_timetag = { 0 }, *Tc_2_timetag = { 0 };  
    int *SW_PVunit = { 0 }, *SW_PV = { 0 }, *SW_Chopper = { 0 };

    cudaMalloc((void**)&VdcBusPU, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&VdcBus, sizeof(double) * 4 * DIM);

    cudaMalloc((void**)&Vcon1A, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Vcon1B, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Vcon1C, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Icon1A, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Icon1B, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Icon1C, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Icon1A1, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Icon1B1, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Icon1C1, sizeof(double) * 4 * DIM);

    cudaMalloc((void**)&Pcon1filt, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Qcon1filt, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&VacRMSfilt, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Ipv, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&VpvBus, sizeof(double) * 4 * DIM);

    cudaMalloc((void**)&VdcPreCh, sizeof(double) * 4 * DIM);

    cudaMalloc((void**)&MPPT_ctrl, sizeof(int) * 4 * DIM);
    cudaMalloc((void**)&Ta_1, sizeof(int) * 4 * DIM);
    cudaMalloc((void**)&Ta_4, sizeof(int) * 4 * DIM);
    cudaMalloc((void**)&Tb_3, sizeof(int) * 4 * DIM);
    cudaMalloc((void**)&Tb_6, sizeof(int) * 4 * DIM);
    cudaMalloc((void**)&Tc_5, sizeof(int) * 4 * DIM);
    cudaMalloc((void**)&Tc_2, sizeof(int) * 4 * DIM);

    cudaMalloc((void**)&MPPT_ctrl_timetag, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Ta_1_timetag, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Ta_4_timetag, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Tb_3_timetag, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Tb_6_timetag, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Tc_5_timetag, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&Tc_2_timetag, sizeof(double) * 4 * DIM);

    cudaMalloc((void**)&SW_PVunit, sizeof(int) * 4 * DIM);
    cudaMalloc((void**)&SW_PV, sizeof(int) * 4 * DIM);
    cudaMalloc((void**)&SW_Chopper, sizeof(int) * 4 * DIM);

    //ATTENTION !!!!!
    double *output_his_LPF_Pcon1 = {0}, double *input_his_LPF_Pcon1 = {0};
    double *output_his_LPF_Qcon1 = {0}, double* input_his_LPF_Qcon1 = {0};
    double *output_his_LPF_VacRMS = {0}, double* input_his_LPF_VacRMS = {0};
    
    cudaMalloc((void**)&output_his_LPF_Pcon1, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&input_his_LPF_Pcon1, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&output_his_LPF_Qcon1, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&input_his_LPF_Qcon1, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&output_his_LPF_VacRMS, sizeof(double) * 4 * DIM);
    cudaMalloc((void**)&input_his_LPF_VacRMS, sizeof(double) * 4 * DIM);

    ////////////////////////Simulation Paramter//////////////////////////////
    double *G_PVunit;
    double *INJ_ABC;

    int *SW_DC;
    int *SW_AC;
    double *Vchar = {0};

    double *d_idc_source;
    double *d_VdcBus;

    double *V_ABC;

    cudaMalloc((void**)&G_PVunit, DIM * 4 * 6 * sizeof(double));
    cudaMalloc((void**)&INJ_ABC, DIM * 4 * 3 * sizeof(double));
   
    cudaMalloc((void**)&SW_DC, DIM * 4 * 4 * sizeof(int));
    //CHECK
    cudaMalloc((void**)&SW_AC, DIM * 4 * 7 * sizeof(int));
    
    //CHECK
    cudaMalloc((void**)&Vchar, DIM * 4 * sizeof(double));

    cudaMalloc((void**)&d_VdcBus, DIM * 4 * sizeof(double));
    cudaMalloc((void**)&d_idc_source, DIM * 4 * sizeof(double));

    cudaMalloc((void**)&V_ABC, DIM * 4 * 3 * sizeof(double));

    //Control static parameters
    /////////////////Start/////////////////
    double* state_his = {0};
    //int* ResetPI1_his = {0}, 
    //int* StartPI1_his = {0};

    cudaMalloc((void**)&state_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&ResetPI1_his, sizeof(int) * DIM * 4);
    //cudaMalloc((void**)&StartPI1_his, sizeof(int) * DIM * 4);

    //LPF_SELF
    double* output_his_LPF_SELF_start = {0};
    double* input_his_LPF_SELF_start = {0};

    cudaMalloc((void**)&output_his_LPF_SELF_start, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&input_his_LPF_SELF_start, sizeof(double) * DIM * 4);

    //SR_Flip_Flop
    int* output_Q_his = {0};
    int* output_Q_not_his = {0};
    int* input_s_his = {0};
    int* input_r_his = {0};

    cudaMalloc((void**)&output_Q_his, sizeof(int) * DIM * 4);
    cudaMalloc((void**)&output_Q_not_his, sizeof(int) * DIM * 4);
    cudaMalloc((void**)&input_s_his, sizeof(int) * DIM * 4);
    cudaMalloc((void**)&input_r_his, sizeof(int) * DIM * 4);

    /////////////////VSC/////////////////
    double* output_his_LPF_VSC = {0};
    double* input_his_LPF_VSC = {0};

    cudaMalloc((void**)&output_his_LPF_VSC, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&input_his_LPF_VSC, sizeof(double) * DIM * 4);
    
    /////////////////PLL and DQ/////////////////
    //double* Vcon1ABpu_his = {0};
    //double* Vcon1BCpu_his = {0}; 
    //double* Vcon1CApu_his = {0};
    //double* Icon1Apu_his = {0};
    //double* Icon1Bpu_his = {0};
    //double* Icon1Cpu_his = {0};

    //cudaMalloc((void**)&Vcon1ABpu_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&Vcon1BCpu_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&Vcon1CApu_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&Icon1Apu_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&Icon1Bpu_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&Icon1Cpu_his, sizeof(double) * DIM * 4);

    //PLL
    PLL3Phase* componentPLL3Phase = {0};
    cudaMalloc((void**)&componentPLL3Phase, DIM * 4 * sizeof(PLL3Phase));

    //LPF
    double* output_his_LPF_PLL = {0};
    double* input_his_LPF_PLL = {0};

    cudaMalloc((void**)&output_his_LPF_PLL, sizeof(double) * DIM * 4 * 4);
    cudaMalloc((void**)&input_his_LPF_PLL, sizeof(double) * DIM * 4 * 4);
    
    /////////////////Voltage Loop/////////////////
    //double* ErrVdcBusPU_his = {0};
    //double* ErrVacPU_his = {0};
    //double* PI_model_output_his1 = {0};
    //double* PI_model_output_his2 = {0};

    //cudaMalloc((void**)&ErrVdcBusPU_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&ErrVacPU_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&PI_model_output_his1, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&PI_model_output_his2, sizeof(double) * DIM * 4);

    //rate_limiter
    double* output_his_rate_limiter_Voltage = {0};
    cudaMalloc((void**)&output_his_rate_limiter_Voltage, sizeof(double) * DIM * 4);

    PI_model* componentPI_model_Voltage = {0};
    cudaMalloc((void**)&componentPI_model_Voltage, DIM * 2 * 4 * sizeof(PI_model));

    /////////////////Current Loop/////////////////
    //double* IdPURef_his = {0};
    //double* IqPURef_his = {0};
    //double* wLeqPU_his = {0};
    //double* Vphcon1PK_his = {0};
    //double* ResetPI1_current_his = {0};

    //cudaMalloc((void**)&IdPURef_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&IqPURef_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&wLeqPU_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&Vphcon1PK_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&ResetPI1_current_his, sizeof(double) * DIM * 4);

    PI_model* componentPI_model_Current = {0};
    cudaMalloc((void**)&componentPI_model_Current, DIM * 2 * 4 * sizeof(PI_model));

    /////////////////Boost/////////////////
    //int* StartPI2_his = {0};
    //int* ResetPI2_his = {0};

    double* output_his_LPF_boost = {0};
    double* input_his_LPF_boost = {0};

    //cudaMalloc((void**)&StartPI2_his, sizeof(int) * DIM * 4);
    //cudaMalloc((void**)&ResetPI2_his, sizeof(int) * DIM * 4);
    cudaMalloc((void**)&output_his_LPF_boost, sizeof(double) * DIM * 4 * 2);
    cudaMalloc((void**)&input_his_LPF_boost, sizeof(double) * DIM * 4 * 2);
    
    MPPT* componentMPPT = {0};
    cudaMalloc((void**)&componentMPPT, DIM * 4 * sizeof(MPPT));

    double* output_his_LPF2_boost = {0};
    double* input_his_LPF2_boost = {0};

    cudaMalloc((void**)&output_his_LPF2_boost, sizeof(double) * DIM * 4);
    cudaMalloc((void**)&input_his_LPF2_boost, sizeof(double) * DIM * 4);

    double* h_output_his_LPF2_boost = {0};
    double* h_input_his_LPF2_boost = {0};

    h_output_his_LPF2_boost = (double*)malloc(sizeof(double) * DIM * 4);
    h_input_his_LPF2_boost = (double*)malloc(sizeof(double) * DIM * 4);

    for(int i = 0; i < DIM * 4; i++)
    {
        h_output_his_LPF2_boost[i] = 2.19759;
        h_input_his_LPF2_boost[i] = 0.0;
    }

    cudaMemcpy(output_his_LPF2_boost, h_output_his_LPF2_boost, sizeof(double) * DIM * 4, cudaMemcpyHostToDevice);
    cudaMemcpy(input_his_LPF2_boost, h_input_his_LPF2_boost, sizeof(double) * DIM * 4, cudaMemcpyHostToDevice);

    PI_model* componentPI_model_boost = {0};
    cudaMalloc((void**)&componentPI_model_boost, DIM * 4 * sizeof(PI_model));

    /////////////////PWM/////////////////
    //double* ModWvA_his = {0};
    //double* ModWvB_his = {0};
    //double* ModWvC_his = {0};

    //cudaMalloc((void**)&ModWvA_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&ModWvB_his, sizeof(double) * DIM * 4);
    //cudaMalloc((void**)&ModWvC_his, sizeof(double) * DIM * 4);

    double* Count_i_Sawtooth_Generator = {0};
    cudaMalloc((void**)&Count_i_Sawtooth_Generator, sizeof(double) * DIM * 4 * 2);

    double* fx_his_Interpolated_Firing_Pulses = {0};
    cudaMalloc((void**)&fx_his_Interpolated_Firing_Pulses, sizeof(double) * DIM * 4 * 7);

    std::ofstream file_PVunit1("data_PVunit1.csv");
    file_PVunit1 << "Time, VpvBus, Ipv, VdcBus, Pcon1filt, Qcon1filt, VacRMSfilt" << endl;

    std::ofstream file_PVunit2("data_PVunit2.csv");
    file_PVunit2 << "Time, VpvBus, Ipv, VdcBus, Pcon1filt, Qcon1filt, VacRMSfilt" << endl;

    std::ofstream file_control("data_control.csv");
    file_control << "Time, SW_PV, SW_Chopper, SW_Char, MPPT_ctrl, Ta_1, Tb_3, Tc_5" << endl;

    std::ofstream file_system("data_system.csv");
    file_system<< "Time, Pstation, Qstation, V1, V2, V3, I1, I2, I3" << endl;

    std::ofstream file_Time("data_Time.csv");
    file_Time<< "Time, control, PVmodel, matrix, system, Total" << endl;
    
    //////////////////////////////////SIMULATION////////////////////////////////
    double *h_simulation_time = { 0 }, *d_simulation_time = { 0 };
    int *h_simulation_num = { 0 }, *d_simulation_num = { 0 };
    double *h_delt = { 0 }, *d_delt = { 0 };

    h_simulation_time = (double*)malloc(sizeof(double));
    h_simulation_num = (int*)malloc(sizeof(int));
    h_delt = (double*)malloc(sizeof(double));

    cudaMalloc((void**)&d_simulation_time, sizeof(double));
    cudaMalloc((void**)&d_simulation_num, sizeof(int));
    cudaMalloc((void**)&d_delt, sizeof(double));

    *h_simulation_time = 0.0;
    *h_simulation_num = 0;
    *h_delt = SIMULATION_DELT;

    cudaMemcpy(d_simulation_time, h_simulation_time, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_simulation_num, h_simulation_num, sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_delt, h_delt, sizeof(double), cudaMemcpyHostToDevice);

    double *h_fault1_time = { 0 }, *h_fault2_time = { 0 };
    double *d_fault1_time = { 0 }, *d_fault2_time = { 0 };
    
    h_fault1_time = (double*)malloc(sizeof(double));
    h_fault2_time = (double*)malloc(sizeof(double));

    cudaMalloc((void**)&d_fault1_time, sizeof(double));
    cudaMalloc((void**)&d_fault2_time, sizeof(double));

    *h_fault1_time = 10.0;
    *h_fault2_time = 4.0;

    cudaMemcpy(d_fault1_time, h_fault1_time, sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_fault2_time, h_fault2_time, sizeof(double), cudaMemcpyHostToDevice);

    //创建cuBLAS句柄
    cublasHandle_t handle;
    cublasCreate(&handle);
    //矩阵参数
    double x = 0.0;
    double y = 1.0;
    double z = -1.0;
    double *beta = &x;
    double *alpha = &y;
    double *alpha1 = &z;
    
    int* info1;
    cudaMalloc(&info1, DIM * 4 * sizeof(int));

    //Total Matrix Calculation
    cusolverDnHandle_t cusolver_handle;
    cusolverDnCreate(&cusolver_handle);
    
    int* d_pivot;
    int* d_info;
    int workspace_size;
	
    cudaMalloc((void**)&d_pivot, sizeof(int) * NODE_NUM);
	cudaMalloc((void**)&d_info, sizeof(int));
    cusolverDnDgetrf_bufferSize(cusolver_handle, NODE_NUM, NODE_NUM, d_Y_total_d, NODE_NUM, &workspace_size);

    double* d_workspace;
    cudaMalloc((void**)&d_workspace, sizeof(double) * workspace_size);

    //data store
    double *h_out_PVunit1_1, *h_out_PVunit1_2, *h_out_PVunit1_3, *h_out_PVunit1_4, *h_out_PVunit1_5, *h_out_PVunit1_6, *h_out_PVunit1_7;
    double *d_out_PVunit1_1, *d_out_PVunit1_2, *d_out_PVunit1_3, *d_out_PVunit1_4, *d_out_PVunit1_5, *d_out_PVunit1_6, *d_out_PVunit1_7;

    double *h_out_PVunit2_1, *h_out_PVunit2_2, *h_out_PVunit2_3, *h_out_PVunit2_4, *h_out_PVunit2_5, *h_out_PVunit2_6, *h_out_PVunit2_7;
    double *d_out_PVunit2_1, *d_out_PVunit2_2, *d_out_PVunit2_3, *d_out_PVunit2_4, *d_out_PVunit2_5, *d_out_PVunit2_6, *d_out_PVunit2_7;

    h_out_PVunit1_1 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	h_out_PVunit1_2 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	h_out_PVunit1_3 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
    h_out_PVunit1_4 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	h_out_PVunit1_5 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	h_out_PVunit1_6 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
    h_out_PVunit1_7 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	
	cudaMalloc((void**)&d_out_PVunit1_1, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_PVunit1_2, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_PVunit1_3, SIMULATION_AMOUNT * sizeof(double));
    cudaMalloc((void**)&d_out_PVunit1_4, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_PVunit1_5, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_PVunit1_6, SIMULATION_AMOUNT * sizeof(double));
    cudaMalloc((void**)&d_out_PVunit1_7, SIMULATION_AMOUNT * sizeof(double));

    h_out_PVunit2_1 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	h_out_PVunit2_2 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	h_out_PVunit2_3 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
    h_out_PVunit2_4 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	h_out_PVunit2_5 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	h_out_PVunit2_6 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
    h_out_PVunit2_7 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	
	cudaMalloc((void**)&d_out_PVunit2_1, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_PVunit2_2, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_PVunit2_3, SIMULATION_AMOUNT * sizeof(double));
    cudaMalloc((void**)&d_out_PVunit2_4, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_PVunit2_5, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_PVunit2_6, SIMULATION_AMOUNT * sizeof(double));
    cudaMalloc((void**)&d_out_PVunit2_7, SIMULATION_AMOUNT * sizeof(double));

    // control data store
    double *h_out_control_1; 
    double *d_out_control_1;
    int *h_out_control_2, *h_out_control_3, *h_out_control_4, *h_out_control_5, *h_out_control_6, *h_out_control_7, *h_out_control_8;
    int *d_out_control_2, *d_out_control_3, *d_out_control_4, *d_out_control_5, *d_out_control_6, *d_out_control_7, *d_out_control_8;

    h_out_control_1 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
    h_out_control_2 = (int*)malloc(SIMULATION_AMOUNT * sizeof(int));
	h_out_control_3 = (int*)malloc(SIMULATION_AMOUNT * sizeof(int));
    h_out_control_4 = (int*)malloc(SIMULATION_AMOUNT * sizeof(int));
	h_out_control_5 = (int*)malloc(SIMULATION_AMOUNT * sizeof(int));
	h_out_control_6 = (int*)malloc(SIMULATION_AMOUNT * sizeof(int));
    h_out_control_7 = (int*)malloc(SIMULATION_AMOUNT * sizeof(int));
    h_out_control_8 = (int*)malloc(SIMULATION_AMOUNT * sizeof(int));
	
	cudaMalloc((void**)&d_out_control_1, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_control_2, SIMULATION_AMOUNT * sizeof(int));
	cudaMalloc((void**)&d_out_control_3, SIMULATION_AMOUNT * sizeof(int));
    cudaMalloc((void**)&d_out_control_4, SIMULATION_AMOUNT * sizeof(int));
	cudaMalloc((void**)&d_out_control_5, SIMULATION_AMOUNT * sizeof(int));
	cudaMalloc((void**)&d_out_control_6, SIMULATION_AMOUNT * sizeof(int));
    cudaMalloc((void**)&d_out_control_7, SIMULATION_AMOUNT * sizeof(int));
    cudaMalloc((void**)&d_out_control_8, SIMULATION_AMOUNT * sizeof(int));

    //sysytem
    double *h_out_system_1, *h_out_system_2, *h_out_system_3, *h_out_system_4, *h_out_system_5, *h_out_system_6, *h_out_system_7, *h_out_system_8, *h_out_system_9;
    double *d_out_system_1, *d_out_system_2, *d_out_system_3, *d_out_system_4, *d_out_system_5, *d_out_system_6, *d_out_system_7, *d_out_system_8, *d_out_system_9;

    h_out_system_1 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	h_out_system_2 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	h_out_system_3 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
    h_out_system_4 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
    h_out_system_5 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	h_out_system_6 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	h_out_system_7 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
    h_out_system_8 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
    h_out_system_9 = (double*)malloc(SIMULATION_AMOUNT * sizeof(double));
	
    cudaMalloc((void**)&d_out_system_1, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_system_2, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_system_3, SIMULATION_AMOUNT * sizeof(double));
    cudaMalloc((void**)&d_out_system_4, SIMULATION_AMOUNT * sizeof(double));
    cudaMalloc((void**)&d_out_system_5, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_system_6, SIMULATION_AMOUNT * sizeof(double));
	cudaMalloc((void**)&d_out_system_7, SIMULATION_AMOUNT * sizeof(double));
    cudaMalloc((void**)&d_out_system_8, SIMULATION_AMOUNT * sizeof(double));
    cudaMalloc((void**)&d_out_system_9, SIMULATION_AMOUNT * sizeof(double));

	//cudaMemcpy(d_out1, h_out1, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyHostToDevice);
	//cudaMemcpy(d_out2, h_out2, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyHostToDevice);
	//cudaMemcpy(d_out3, h_out3, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyHostToDevice);
    //cudaMemcpy(d_out4, h_out4, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyHostToDevice);
	//cudaMemcpy(d_out5, h_out5, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyHostToDevice);
	//cudaMemcpy(d_out6, h_out6, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyHostToDevice);
    //cudaMemcpy(d_out7, h_out7, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyHostToDevice);

    //计时
    cudaEvent_t time_total1, time_total2;
    float elapsedTime_total;
    cudaEvent_t time_control1, time_control2, time_pvmodel1, time_pvmodel2, time_matrix1, time_matrix2, time_system1, time_system2;
    float elapsedTime_control, elapsedTime_pvmodel, elapsedTime_matrix, elapsedTime_system;

    //system
    //cudaEvent_t time_system_model1, time_system_model2, time_system_matrix1, time_system_matrix2;
    //cudaEvent_t time_system_YJtotal1, time_system_YJtotal2, time_system_Vtotal1, time_system_Vtotal2, time_system_output1, time_system_output2;
    //float elapsedTime_system_model, elapsedTime_system_matrix, elapsedTime_system_YJtotal, elapsedTime_system_Vtotal, elapsedTime_system_output;

    float *TimeStore_total;
    float *TimeStore_control = { 0 }, *TimeStore_pvmodel = { 0 }, *TimeStore_matrix = { 0 }, *TimeStore_system = { 0 };
    //float *TimeStore_system_model = { 0 }, *TimeStore_system_matrix = { 0 };
    //float *TimeStore_system_YJtotal = { 0 }, *TimeStore_system_Vtotal = { 0 }, *TimeStore_system_output = { 0 };
    
    TimeStore_total = (float*)malloc(SIMULATION_AMOUNT * sizeof(float));
    TimeStore_control = (float*)malloc(SIMULATION_AMOUNT * sizeof(float));
    TimeStore_pvmodel = (float*)malloc(SIMULATION_AMOUNT * sizeof(float));
    TimeStore_matrix = (float*)malloc(SIMULATION_AMOUNT * sizeof(float));
    TimeStore_system = (float*)malloc(SIMULATION_AMOUNT * sizeof(float));

    //TimeStore_system_model = (float*)malloc(SIMULATION_AMOUNT * sizeof(float));
    //TimeStore_system_matrix = (float*)malloc(SIMULATION_AMOUNT * sizeof(float));
    //TimeStore_system_YJtotal = (float*)malloc(SIMULATION_AMOUNT * sizeof(float));
    //TimeStore_system_Vtotal = (float*)malloc(SIMULATION_AMOUNT * sizeof(float));
    //TimeStore_system_output = (float*)malloc(SIMULATION_AMOUNT * sizeof(float));

    cudaEventCreate(&time_total1);
    cudaEventCreate(&time_total2);
    cudaEventCreate(&time_control1);
    cudaEventCreate(&time_control2);
    cudaEventCreate(&time_pvmodel1);
    cudaEventCreate(&time_pvmodel2);
    cudaEventCreate(&time_matrix1);
    cudaEventCreate(&time_matrix2);
    cudaEventCreate(&time_system1);
    cudaEventCreate(&time_system2);

    //cudaEventCreate(&time_system_model1);
    //cudaEventCreate(&time_system_model2);
    //cudaEventCreate(&time_system_matrix1);
    //cudaEventCreate(&time_system_matrix2);
    //cudaEventCreate(&time_system_YJtotal1);
    //cudaEventCreate(&time_system_YJtotal2);
    //cudaEventCreate(&time_system_Vtotal1);
    //cudaEventCreate(&time_system_Vtotal2);
    //cudaEventCreate(&time_system_output1);
    //cudaEventCreate(&time_system_output2);

    //流处理
    cudaStream_t stream_pvmodel1, stream_pvmodel2;
    cudaStreamCreate(&stream_pvmodel1);
    cudaStreamCreate(&stream_pvmodel2);

    cudaStream_t stream_batch1, stream_batch2, stream_batch3, stream_batch4;
    cudaStreamCreate(&stream_batch1);
    cudaStreamCreate(&stream_batch2);
    cudaStreamCreate(&stream_batch3);
    cudaStreamCreate(&stream_batch4);

    cudaStream_t stream_matrix1, stream_matrix2, stream_matrix3, stream_matrix4, stream_matrix5, stream_matrix6;
    cudaStreamCreate(&stream_matrix1);
    cudaStreamCreate(&stream_matrix2);
    cudaStreamCreate(&stream_matrix3);
    cudaStreamCreate(&stream_matrix4);
    cudaStreamCreate(&stream_matrix5);
    cudaStreamCreate(&stream_matrix6);

    while(*h_simulation_time < TOTAL_TIME)
    {  
        cudaEventRecord(time_total1, 0);
        cudaEventRecord(time_control1, 0);
        
        /////////////////Control Calculation/////////////////
        //修改
        //Dblk[i] = 0; //0.08s之前是0之后是1；比较器time和0.08s
        //SW_Char[i] = 0; //0.05s之前是0之后是1;比较器time和0.05s
        Control_Init << < 1, DIM * 4 >> > (d_simulation_time, d_Dblk, d_SW_Char);
        cudaDeviceSynchronize();

        Control_PVmodel_kernel << <1, DIM * 4>> > (d_Dblk, d_comparator_1_B, d_comparator_2_A, d_VLLrated, d_fsys, d_Lcon1, d_Vtr2Pri, d_Vtr2Sec, d_PoRated, 
        Vcon1A, Vcon1B, Vcon1C, Icon1A1, Icon1B1, Icon1C1, VdcBus, d_VdcRated, Pcon1filt, Qcon1filt, VacRMSfilt, Ipv, VpvBus, 
        d_VdcCon, d_VpvRated, d_signal_A, d_FrzPI, d_VdcPURef, d_VacPURef, d_QPURef, d_VacOrQ, d_IdMax, d_IdMin, d_Park0, d_Har3ctrl, d_fsw1, d_fsw2, VdcPreCh, 
        MPPT_ctrl, MPPT_ctrl_timetag, Ta_1, Ta_1_timetag, Ta_4, Ta_4_timetag, Tb_3, Tb_3_timetag, Tb_6, Tb_6_timetag, Tc_5, Tc_5_timetag, Tc_2, Tc_2_timetag, 
        d_G_VSC, d_T_VSC, d_Strtup, d_Samp, d_Vstep, d_Sd, d_G1_Boost, d_T1_Boost, d_G2_Boost, d_T2_Boost, d_G3_Boost, d_T3_Boost, d_signal, d_Max_Boost, d_Min_Boost, d_Kp_Boost, d_Ti_Boost, d_Kp, d_Kp2, d_ULIM, d_LLIM, d_MAXOUT, d_MINOUT, d_MAXOUT2, d_MINOUT2, d_PI_T1, d_PI_T2, 
        d_Mmax_current, d_Mmin_current, d_MaxLimit1_current, d_MinLimit1_current, d_MaxLimit2_current, d_MinLimit2_current, d_Kp1_current, d_Kp2_current, d_Ti1_current, d_Ti2_current, 
        d_Kp_PLL, d_Ti_PLL, d_Base_Vpll, d_fbase, d_fmaxPU, d_fminPU, d_Kp_realpole, d_Ti_realpole_L, d_Ti_realpole_R, d_Ctrl, d_max1_pwm, d_min1_pwm, d_max2_pwm, d_min2_pwm, 
        d_G_start, d_T_start, d_comparator1_H, d_comparator1_L, d_comparator2_H, d_comparator2_L, SW_PVunit, SW_PV, SW_Chopper, d_delt, d_simulation_time, d_theta1, d_theta2,
        state_his, output_his_LPF_SELF_start, input_his_LPF_SELF_start, output_Q_his, output_Q_not_his, input_s_his, input_r_his, output_his_LPF_VSC, input_his_LPF_VSC,
        output_his_LPF_PLL, input_his_LPF_PLL, componentPLL3Phase,
        output_his_rate_limiter_Voltage, componentPI_model_Voltage,
        componentPI_model_Current,
        output_his_LPF_boost, input_his_LPF_boost, output_his_LPF2_boost, input_his_LPF2_boost, componentMPPT, componentPI_model_boost,
        Count_i_Sawtooth_Generator, fx_his_Interpolated_Firing_Pulses);
        cudaDeviceSynchronize();

        SW_AC_DC << < 1, DIM * 4 >> > (SW_PV, SW_Chopper, d_SW_Char, MPPT_ctrl, Ta_1, Tb_3, Tc_5, Ta_4, Tb_6, Tc_2, SW_PVunit, VdcPreCh,
        SW_DC, SW_AC, Vchar);
        cudaDeviceSynchronize();

        cudaEventRecord(time_control2, 0);
        cudaEventSynchronize(time_control2);

        cudaEventRecord(time_pvmodel1, 0);

        if(*h_simulation_time == 0)
        {
            pvmodel_pre_dc << < 1, DIM * 4>> > (d_delt, d_PVarray_pars, d_Boost_pars, d_CPB_pars, d_VSC_pars, 
                Bre_PV, Bre_Chop, Bre_Char, Boost_IGBT, Boost_Diode, C_PV, C_Boost, C_VSC, L_Boost,
                Npar, Nser, Isc1_PV, Voc1_PV, par_PV_c1, par_PV_b1, V_PV_V_his, I_PV_A, I_PV_kA, STATUS_his,
                d_Yg2_d, d_J_branch2_d);

            cudaDeviceSynchronize();
            
            pvmodel_pre_ac << < 1, DIM * 4 >> > (d_delt, d_VSC_pars, d_LC_pars, d_PVunit_Bre_pars, d_Trans_pars,
                Bre_PVunit, VSC_average, C_filt, L_filt, 
                G_filt_R1, G_filt_R2, d_G_trans, d_parB_trans, d_J_trans_A_PVunit_d, d_J_trans_B_PVunit_d, d_J_trans_C_PVunit_d,
                d_Yg1_d, d_J_branch1_d);

            cudaDeviceSynchronize();    

            //cudaStreamSynchronize(stream_pvmodel1);
            //cudaStreamSynchronize(stream_pvmodel2);
        }
        else
        {
            //反解
            //Jin - Y21 * Vex
            cublasDgemmBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, NODE_NUM_AC -3, 1, 3, alpha, d_Y21_d, NODE_NUM_AC - 3, d_Vex_d, 3, beta, d_YV_d, NODE_NUM_AC - 3, DIM * 4);
            cudaDeviceSynchronize();

            JYV_Calculation << < 1, DIM * 4>> > (d_Jin_d, d_YV_d, d_JYV_d);
            cudaDeviceSynchronize();

            //Vin
            cublasDgemmBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, NODE_NUM_AC - 3, 1, NODE_NUM_AC - 3, alpha, d_Y22_inv_d, NODE_NUM_AC - 3, d_JYV_d, NODE_NUM_AC - 3, beta, d_Vin_d, NODE_NUM_AC - 3, DIM * 4);
            cudaDeviceSynchronize();
            
            pvmodel_dc << < 1, DIM * 4 >> > (SW_DC, Vchar, d_idc_source, d_V_DC_d,
                Bre_PV, Bre_Chop, Bre_Char, Boost_IGBT, Boost_Diode, C_PV, C_Boost, C_VSC, L_Boost,
                Npar, Nser, Isc1_PV, Voc1_PV, par_PV_c1, par_PV_b1, V_PV_V_his, V_PV_V, V_PV_kV, I_PV_A, I_PV_kA, STATUS_his,
                d_Yg2_d, d_J_branch2_d);

            cudaDeviceSynchronize();

            pvmodel_ac << < 1, DIM * 4 >> > (SW_AC, d_Dblk, VdcBus, d_idc_source,
                d_Vex_d, d_Vin_d, I_port_A_PVunit, I_port_B_PVunit, I_port_C_PVunit,
                Bre_PVunit, VSC_average, C_filt, L_filt, 
                G_filt_R1, G_filt_R2, d_G_trans, d_parB_trans, d_J_trans_A_PVunit_d, d_J_trans_B_PVunit_d, d_J_trans_C_PVunit_d,
                d_Yg1_d, d_J_branch1_d);

            cudaDeviceSynchronize();
            
            //cudaStreamSynchronize(stream_pvmodel1);
            //cudaStreamSynchronize(stream_pvmodel2);         
        }

        cudaEventRecord(time_pvmodel2, 0);
        cudaEventSynchronize(time_pvmodel2);

        cudaEventRecord(time_matrix1, 0);

        //Matrix Calculation
        //Gernerate Y and J
        //AC side
        //Y1 = A1 * Yg1 * A1.transpose(); alpha = 1
        //cublasSetStream(handle, stream_batch1);
        cublasDgemmBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, NODE_NUM_AC, BRANCH_NUM_AC, BRANCH_NUM_AC, alpha, d_A1_d, NODE_NUM_AC, d_Yg1_d, BRANCH_NUM_AC, beta, d_AYg1_d, NODE_NUM_AC, DIM * 4);
        cudaDeviceSynchronize();
        cublasDgemmBatched(handle, CUBLAS_OP_N, CUBLAS_OP_T, NODE_NUM_AC, NODE_NUM_AC, BRANCH_NUM_AC, alpha, d_AYg1_d, NODE_NUM_AC, d_A1_d, NODE_NUM_AC, beta, d_Y1_d, NODE_NUM_AC, DIM * 4);
        cudaDeviceSynchronize();
        
        //J1 = - A1 * J_branch1; alpha1 = -1
        //cublasSetStream(handle, stream_batch2);
        cublasDgemmBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, NODE_NUM_AC, 1, BRANCH_NUM_AC, alpha1, d_A1_d, NODE_NUM_AC, d_J_branch1_d, BRANCH_NUM_AC, beta, d_J1_d, NODE_NUM_AC, DIM * 4);
        cudaDeviceSynchronize();

        //DC side
        //Y2 = A2 * Yg2 * A2.transpose(); alpha = 1
        //cublasSetStream(handle, stream_batch3);
        cublasDgemmBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, NODE_NUM_DC, BRANCH_NUM_DC, BRANCH_NUM_DC, alpha, d_A2_d, NODE_NUM_DC, d_Yg2_d, BRANCH_NUM_DC, beta, d_AYg2_d, NODE_NUM_DC, DIM * 4);
        cudaDeviceSynchronize();
        cublasDgemmBatched(handle, CUBLAS_OP_N, CUBLAS_OP_T, NODE_NUM_DC, NODE_NUM_DC, BRANCH_NUM_DC, alpha, d_AYg2_d, NODE_NUM_DC, d_A2_d, NODE_NUM_DC, beta, d_Y2_d, NODE_NUM_DC, DIM * 4);
        cudaDeviceSynchronize();

        //J2 = - A2 * J_branch2; alpha1 = -1
        //cublasSetStream(handle, stream_batch4);
        cublasDgemmBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, NODE_NUM_DC, 1, BRANCH_NUM_DC, alpha1, d_A2_d, NODE_NUM_DC, d_J_branch2_d, BRANCH_NUM_DC, beta, d_J2_d, NODE_NUM_DC, DIM * 4);
        cudaDeviceSynchronize();
        
        //cudaStreamSynchronize(stream_batch1);
        //cudaStreamSynchronize(stream_batch2);
        //cudaStreamSynchronize(stream_batch3);
        //cudaStreamSynchronize(stream_batch4);

        //矩阵分块
        int threadsPerBlock = 256;
        int blocksPerGrid = DIM * 4;  // 每个网格处理一个矩阵
        MatrixPartitioning_Y11<<<blocksPerGrid, threadsPerBlock>>>(d_Y1_d, d_Y11_d);
        MatrixPartitioning_Y12<<<blocksPerGrid, threadsPerBlock>>>(d_Y1_d, d_Y12_d);
        MatrixPartitioning_Y21<<<blocksPerGrid, threadsPerBlock>>>(d_Y1_d, d_Y21_d);
        MatrixPartitioning_Y22<<<blocksPerGrid, threadsPerBlock>>>(d_Y1_d, d_Y22_d);
        MatrixPartitioning_Jex<<<blocksPerGrid, threadsPerBlock>>>(d_J1_d, d_Jex_d);
        MatrixPartitioning_Jin<<<blocksPerGrid, threadsPerBlock>>>(d_J1_d, d_Jin_d);

        cudaDeviceSynchronize();
        
        //cudaStreamSynchronize(stream_matrix1);
        //cudaStreamSynchronize(stream_matrix2);
        //cudaStreamSynchronize(stream_matrix3);
        //cudaStreamSynchronize(stream_matrix4);
        //cudaStreamSynchronize(stream_matrix5);
        //cudaStreamSynchronize(stream_matrix6);

        //Y22_inv
        //int* info1;
        //cudaMalloc(&info1, DIM * sizeof(int));
        cublasDmatinvBatched(handle, (NODE_NUM_AC - 3), d_Y22_d, (NODE_NUM_AC - 3), d_Y22_inv_d, (NODE_NUM_AC - 3), info1, DIM * 4);
        cudaDeviceSynchronize();

        //Y12 * Y22_inv
        cublasDgemmBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, 3, NODE_NUM_AC - 3, NODE_NUM_AC - 3, alpha, d_Y12_d, 3, d_Y22_inv_d, NODE_NUM_AC - 3, beta, d_YY_d, 3, DIM * 4);
        cudaDeviceSynchronize();

        //Y12 * Y22_inv * Y21
        //cublasSetStream(handle, stream_batch1);
        cublasDgemmBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, 3, 3, NODE_NUM_AC - 3, alpha, d_YY_d, 3, d_Y21_d, NODE_NUM_AC - 3, beta, d_YYY_d, 3, DIM * 4);
        cudaDeviceSynchronize();
        
        //Y12 * Y22_inv * Jin
        //cublasSetStream(handle, stream_batch2);
        cublasDgemmBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, 3, 1, NODE_NUM_AC - 3, alpha, d_YY_d, 3, d_Jin_d, NODE_NUM_AC - 3, beta, d_YYJ_d, 3, DIM * 4);
        cudaDeviceSynchronize();
        
        //DC side
        //V_DC = Y2.inverse() * J2
        //cublasSetStream(handle, stream_batch3);
        cublasDmatinvBatched(handle, NODE_NUM_DC, d_Y2_d, NODE_NUM_DC, d_Y2_inv_d, NODE_NUM_DC, info1, DIM * 4);
        cudaDeviceSynchronize();
        cublasDgemmBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, NODE_NUM_DC, 1, NODE_NUM_DC, alpha, d_Y2_inv_d, NODE_NUM_DC, d_J2_d, NODE_NUM_DC, beta, d_V_DC_d, NODE_NUM_DC, DIM * 4);
        cudaDeviceSynchronize();
        
        //cudaStreamSynchronize(stream_batch1);
        //cudaStreamSynchronize(stream_batch2);
        //cudaStreamSynchronize(stream_batch3);

        //Yeq & Jeq Calculation
        Yeq_Jeq_Calculation << < 1, DIM * 4 >> > (d_Y11_d, d_YYY_d, d_Yeq_d, d_Jex_d, d_YYJ_d, d_Jeq_d, G_PVunit, INJ_ABC);
        cudaDeviceSynchronize();

        cudaEventRecord(time_matrix2, 0);
        cudaEventSynchronize(time_matrix2);

        cudaEventRecord(time_system1, 0);

        //cudaEventRecord(time_system_model1, 0);

        //System Simulation
        if(*h_simulation_time == 0)
        {
            system_pre <<<1, 1>>> (d_delt, d_Source_pars, d_Fault_pars, d_PImodel_1_pars, d_Trans_sys_pars, d_PImodel_2_pars,
                L_source, Fault_1, Fault_2, PImodel_1, PImodel_2, d_G_trans_sys, d_parB_trans_sys);

            cudaDeviceSynchronize();
        }
        
        system_emt <<<1, 1>>>(d_Source_pars, d_simulation_time, d_fault1_time, d_fault2_time, L_source, Fault_1, Fault_2, PImodel_1, PImodel_2,
            d_V_total_d, d_Yg_CAL_d, d_J_branch_CAL_d, G_PVunit, INJ_ABC,
            d_G_trans_sys, d_parB_trans_sys, d_J_trans_A_sys_d, d_J_trans_B_sys_d, d_J_trans_C_sys_d,
            d_Ih_mutual1_1, d_Ih_mutual2_1, d_Ih_self1_1, d_Ih_self2_1, d_Ih_ABC_1,
            d_Ih_mutual1_2, d_Ih_mutual2_2, d_Ih_self1_2, d_Ih_self2_2, d_Ih_ABC_2,
            d_Pstation, output_his_LPF_Pstation, input_his_LPF_Pstation, d_Qstation, output_his_LPF_Qstation, input_his_LPF_Qstation, d_Vstation, d_Istation);

        cudaDeviceSynchronize();

        //cudaEventRecord(time_system_model2, 0);
        //cudaEventSynchronize(time_system_model2);

        //cudaEventRecord(time_system_matrix1, 0);
        
	    //Y_CAL = A_CAL * Yg_CAL * A_CAL.transpose();
        cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, NODE_NUM, BRANCH_NUM, BRANCH_NUM, alpha, d_A_CAL_d, NODE_NUM, d_Yg_CAL_d, BRANCH_NUM, beta, d_AYg_CAL_d, NODE_NUM);
        cudaDeviceSynchronize();
        cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, NODE_NUM, NODE_NUM, BRANCH_NUM, alpha, d_AYg_CAL_d, NODE_NUM, d_A_CAL_d, NODE_NUM, beta, d_Y_CAL_d, NODE_NUM);
        cudaDeviceSynchronize();

        //J_CAL = - A_CAL * J_branch_CAL;
        cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, NODE_NUM, 1, BRANCH_NUM, alpha1, d_A_CAL_d, NODE_NUM, d_J_branch_CAL_d, BRANCH_NUM, beta, d_J_CAL_d, NODE_NUM);
        cudaDeviceSynchronize();

        //cudaEventRecord(time_system_matrix2, 0);
        //cudaEventSynchronize(time_system_matrix2);

        //cudaEventRecord(time_system_YJtotal1, 0);
        
        //Y_total Calculation
        Y_total_Calculation <<<1, 1>>>(d_Y_total_d, d_Y_total_inv_d, d_Y_CAL_d, PImodel_1, PImodel_2);
        cudaDeviceSynchronize();

        //J_total Calculation
        J_total_Calculation <<<1, 1>>>(d_J_total_d, d_J_CAL_d,
            d_Ih_mutual1_1, d_Ih_mutual2_1, d_Ih_self1_1, d_Ih_self2_1, d_Ih_ABC_1,
            d_Ih_mutual1_2, d_Ih_mutual2_2, d_Ih_self1_2, d_Ih_self2_2, d_Ih_ABC_2);
        cudaDeviceSynchronize();

        //cudaEventRecord(time_system_YJtotal2, 0);
        //cudaEventSynchronize(time_system_YJtotal2);

        //cudaEventRecord(time_system_Vtotal1, 0);
    
        //V_total = Y_total.inverse() * J_total;
        //LU //CHECK!!!!
        cusolverDnDgetrf(cusolver_handle, NODE_NUM, NODE_NUM, d_Y_total_d, NODE_NUM, d_workspace, d_pivot, d_info);
        cudaDeviceSynchronize();

        cusolverDnDgetrs(cusolver_handle, CUBLAS_OP_N, NODE_NUM, NODE_NUM, d_Y_total_d, NODE_NUM, d_pivot, d_Y_total_inv_d, NODE_NUM, d_info);
        cudaDeviceSynchronize();

        //size_t shared_mem_size = 2 * NODE_NUM * NODE_NUM * sizeof(double); // 共享内存大小 
        //invertMatrixKernel<<<1, 32, shared_mem_size>>>(d_Y_total_d, d_Y_total_inv_d, NODE_NUM); 
        //cudaDeviceSynchronize;

        cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, NODE_NUM, 1, NODE_NUM, alpha, d_Y_total_inv_d, NODE_NUM, d_J_total_d, NODE_NUM, beta, d_V_total_d, NODE_NUM);
        cudaDeviceSynchronize();

        //cudaEventRecord(time_system_Vtotal2, 0);
        //cudaEventSynchronize(time_system_Vtotal2);

        //cudaEventRecord(time_system_output1, 0);

        //PVunit output
        //System
        Output_System <<< 1 , DIM >>>(d_V_total_d, V_ABC, d_Vex_d);
        cudaDeviceSynchronize();
        
        //DCside
        Output_DCside <<< 1, DIM * 4>>>(V_PV_kV, I_PV_kA, d_V_DC_d, VpvBus, Ipv, VdcBus);
        cudaDeviceSynchronize();
        
        //ACside
        Output_ACside <<< 1, DIM * 4>>>(d_Vin_d, I_port_A_PVunit, I_port_B_PVunit, I_port_C_PVunit,
        Vcon1A, Vcon1B, Vcon1C, Icon1A, Icon1B, Icon1C, Icon1A1, Icon1B1, Icon1C1,
        output_his_LPF_Pcon1, input_his_LPF_Pcon1, output_his_LPF_Qcon1, input_his_LPF_Qcon1, output_his_LPF_VacRMS, input_his_LPF_VacRMS, 
        Pcon1filt, Qcon1filt, VacRMSfilt);

        cudaDeviceSynchronize();

        Output_print <<< 1, 1 >>> (d_simulation_time, d_simulation_num, VpvBus, Ipv, VdcBus, Pcon1filt, Qcon1filt, VacRMSfilt, SW_AC, SW_DC, d_Pstation, d_Qstation, d_Vstation, d_Istation,
        d_out_PVunit1_1, d_out_PVunit1_2, d_out_PVunit1_3, d_out_PVunit1_4, d_out_PVunit1_5, d_out_PVunit1_6, d_out_PVunit1_7,
        d_out_PVunit2_1, d_out_PVunit2_2, d_out_PVunit2_3, d_out_PVunit2_4, d_out_PVunit2_5, d_out_PVunit2_6, d_out_PVunit2_7,
        d_out_control_1, d_out_control_2, d_out_control_3, d_out_control_4, d_out_control_5, d_out_control_6, d_out_control_7, d_out_control_8,
        d_out_system_1, d_out_system_2, d_out_system_3, d_out_system_4, d_out_system_5, d_out_system_6, d_out_system_7, d_out_system_8, d_out_system_9);
        
        cudaDeviceSynchronize();

        //cudaEventRecord(time_system_output2, 0);
        //cudaEventSynchronize(time_system_output2);

        cudaEventRecord(time_system2, 0);
        cudaEventSynchronize(time_system2);

        cudaEventRecord(time_total2, 0);
        cudaEventSynchronize(time_total2);

        cudaEventElapsedTime(&elapsedTime_total, time_total1, time_total2);
        cudaEventElapsedTime(&elapsedTime_control, time_control1, time_control2);
        cudaEventElapsedTime(&elapsedTime_pvmodel, time_pvmodel1, time_pvmodel2);
        cudaEventElapsedTime(&elapsedTime_matrix, time_matrix1, time_matrix2);
        cudaEventElapsedTime(&elapsedTime_system, time_system1, time_system2);

        //cudaEventElapsedTime(&elapsedTime_system_model, time_system_model1, time_system_model2);
        //cudaEventElapsedTime(&elapsedTime_system_matrix, time_system_matrix1, time_system_matrix2);
        //cudaEventElapsedTime(&elapsedTime_system_YJtotal, time_system_YJtotal1, time_system_YJtotal2);
        //cudaEventElapsedTime(&elapsedTime_system_Vtotal, time_system_Vtotal1, time_system_Vtotal2);
        //cudaEventElapsedTime(&elapsedTime_system_output, time_system_output1, time_system_output2);

        TimeStore_total[*h_simulation_num] = elapsedTime_total;
        TimeStore_control[*h_simulation_num] = elapsedTime_control;
        TimeStore_pvmodel[*h_simulation_num] = elapsedTime_pvmodel;
        TimeStore_matrix[*h_simulation_num] = elapsedTime_matrix;
        TimeStore_system[*h_simulation_num] = elapsedTime_system;

        //TimeStore_system_model[*h_simulation_num] = elapsedTime_system_model;
        //TimeStore_system_matrix[*h_simulation_num] = elapsedTime_system_matrix;
        //TimeStore_system_YJtotal[*h_simulation_num] = elapsedTime_system_YJtotal;
        //TimeStore_system_Vtotal[*h_simulation_num] = elapsedTime_system_Vtotal;
        //TimeStore_system_output[*h_simulation_num] = elapsedTime_system_output;

        //cout << "Time:" << TimeStore[*h_simulation_num] << endl;

        *h_simulation_num += 1;
        *h_simulation_time += SIMULATION_DELT;
    }

    cudaMemcpy(h_out_PVunit1_1, d_out_PVunit1_1, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_PVunit1_2, d_out_PVunit1_2, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_PVunit1_3, d_out_PVunit1_3, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
    cudaMemcpy(h_out_PVunit1_4, d_out_PVunit1_4, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_PVunit1_5, d_out_PVunit1_5, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_PVunit1_6, d_out_PVunit1_6, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
    cudaMemcpy(h_out_PVunit1_7, d_out_PVunit1_7, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);

    cudaMemcpy(h_out_PVunit2_1, d_out_PVunit2_1, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_PVunit2_2, d_out_PVunit2_2, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_PVunit2_3, d_out_PVunit2_3, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
    cudaMemcpy(h_out_PVunit2_4, d_out_PVunit2_4, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_PVunit2_5, d_out_PVunit2_5, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_PVunit2_6, d_out_PVunit2_6, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
    cudaMemcpy(h_out_PVunit2_7, d_out_PVunit2_7, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);

    cudaMemcpy(h_out_control_1, d_out_control_1, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_control_2, d_out_control_2, SIMULATION_AMOUNT * sizeof(int), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_control_3, d_out_control_3, SIMULATION_AMOUNT * sizeof(int), cudaMemcpyDeviceToHost);
    cudaMemcpy(h_out_control_4, d_out_control_4, SIMULATION_AMOUNT * sizeof(int), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_control_5, d_out_control_5, SIMULATION_AMOUNT * sizeof(int), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_control_6, d_out_control_6, SIMULATION_AMOUNT * sizeof(int), cudaMemcpyDeviceToHost);
    cudaMemcpy(h_out_control_7, d_out_control_7, SIMULATION_AMOUNT * sizeof(int), cudaMemcpyDeviceToHost);
    cudaMemcpy(h_out_control_8, d_out_control_8, SIMULATION_AMOUNT * sizeof(int), cudaMemcpyDeviceToHost);

    cudaMemcpy(h_out_system_1, d_out_system_1, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_system_2, d_out_system_2, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_system_3, d_out_system_3, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
    cudaMemcpy(h_out_system_4, d_out_system_4, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
    cudaMemcpy(h_out_system_5, d_out_system_5, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_system_6, d_out_system_6, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(h_out_system_7, d_out_system_7, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
    cudaMemcpy(h_out_system_8, d_out_system_8, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
    cudaMemcpy(h_out_system_9, d_out_system_9, SIMULATION_AMOUNT * sizeof(double), cudaMemcpyDeviceToHost);
    
    for(int i = 0; i < SIMULATION_AMOUNT; i++)
    {
        file_PVunit1 << h_out_PVunit1_1[i] << ',' << h_out_PVunit1_2[i] << ',' << h_out_PVunit1_3[i] << ',' << h_out_PVunit1_4[i] << ',';
        file_PVunit1 << h_out_PVunit1_5[i] << ',' << h_out_PVunit1_6[i] << ',' << h_out_PVunit1_7[i] << endl;

        file_PVunit2 << h_out_PVunit2_1[i] << ',' << h_out_PVunit2_2[i] << ',' << h_out_PVunit2_3[i] << ',' << h_out_PVunit2_4[i] << ',';
        file_PVunit2 << h_out_PVunit2_5[i] << ',' << h_out_PVunit2_6[i] << ',' << h_out_PVunit2_7[i] << endl;

        file_control << h_out_control_1[i] << ',' ;
        file_control << h_out_control_2[i] << ',' << h_out_control_3[i] << ',' << h_out_control_4[i] << ',' << h_out_control_5[i] << ',';
        file_control << h_out_control_6[i] << ',' << h_out_control_7[i] << ',' << h_out_control_8[i] << endl;

        file_system << h_out_system_1[i] << ',' << h_out_system_2[i] << ',' << h_out_system_3[i] << ',';
        file_system << h_out_system_4[i] << ',' << h_out_system_5[i] << ',' << h_out_system_6[i] << ',';
        file_system << h_out_system_7[i] << ',' << h_out_system_8[i] << ',' << h_out_system_9[i] << endl;

        file_Time << h_out_control_1[i] << ',' << TimeStore_control[i] << ',' << TimeStore_pvmodel[i] << ',';
        file_Time << TimeStore_matrix[i] << ',' << TimeStore_system[i] << ',' << TimeStore_total[i] << endl;
    }

    file_PVunit1.close();
    //file3.close();
    //file4.close();

    cudaStreamDestroy(stream_pvmodel1);
    cudaStreamDestroy(stream_pvmodel2);
    
    cudaStreamDestroy(stream_matrix1);
    cudaStreamDestroy(stream_matrix2);
    cudaStreamDestroy(stream_matrix3);
    cudaStreamDestroy(stream_matrix4);
    cudaStreamDestroy(stream_matrix5);
    cudaStreamDestroy(stream_matrix6);

    cudaStreamDestroy(stream_batch1);
    cudaStreamDestroy(stream_batch2);
    cudaStreamDestroy(stream_batch3);
    cudaStreamDestroy(stream_batch4);
    
    
    return 0;

}